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ABSTRACT 


Hie  axisymmetric  heat,  equation,  resulting  from  a  point-souce  of  heat  applied 
lo  a  metal  block,  is  solved  numerically;  both  iterative  and  multilevel  solutions  are 
computed  in  order  to  compare  the  two  processes.  The  continuum  problem  is  dis- 
cretized  in  two  stages:  finite  differences  are  used  to  discretize  the  time  derivatives, 
resulting  in  a  fully  implicit  backward  time-stepping  scheme,  and  the  Finite  Volume 
Element  (FVE)  method  is  used  to  discretize  the  spatial  derivatives.  The  application 
of  tin1  FV'E  method  to  a  problem  in  cylindrical  coordinates  is  new.  and  results  in 
stencils  which  are  analyzed  extensively.  Several  iteration  schemes  are  considered,  in- 
<  hiding  both  Jacobi  and  Gauss-Seidel;  a  thorough  analysis  of  these  schemes  is  done, 
using  both  the  spectral  radii  of  the  iteration  matrices  and  local  mode  analysis.  Us- 
ing this  discretization,  a  Gauss-Seidel  relaxation  scheme  is  used  to  solve  the  heat 
equation  iteratively.  A  multilevel  solution  process  is  then  constructed,  including  the 
development  of  intergrid  transfer  and  coarse  grid  operators.  Local  mode  analysis  is 
performed  on  the  components  of  the  amplification  matrix,  resulting  in  the  two-level 
convergence  factors  for  various  combinations  of  the  operators.  The  multilevel  solu- 
t  ion  process  is  implemented  by  using  multigrid  V-cycles;  the  iterative  and  multilevel 
results  are  compared  and  discussed  in  detail.  The  computational  savings  resulting 
from  the  multilevel  process  are  then  discussed. 
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I.         INTRODUCTION 


A  topic  oi  genera]  interest  is  the  numerical  solution  of  partial  differential  equa- 
tions (PDEs),  such  as  the  Navier-Stokes  equation  (Equation  II. 1).  Many  aspects  of 
these  types  oi  problems  present  interesting  challenges  to  constructing  numerical  so- 
lutions, such  as  non-linearities,  high-speed  flows  that  cause  numerical  peculiarities, 
unusual  boundary  conditions,  uncommon  geometries,  and  so  on.  Developing  pro- 
cesses that  solve  these  types  of  problems  is  difficult,  and  often  results  in  complicated 
solutions  schemes  which  are  costly  to  implement.  A  common  goal  is  to  attempt  to 
streamline  existing  solution  processes,  or  develop  new  ones,  in  order  to  reduce  com- 
putational complexity  and  cost. 

At  issue  is  how  to  transform  a  continuum  problem  into  a  discrete  one  (dis- 
cretization), and  how  to  construct  a  numerical  process  that  will  solve  the  discrete 
problem.  The  finite  volume  element  (FVE)  method  (see  [Ref.  1])  has  proven  to  be 
a  useful  tool  in  developing  discretizations,  particularly  for  problems  that  require  the 
enforcement  of  conservation  laws.  One  of  the  more  successful  methods  for  streamlin- 
ing the  solution  of  PDEs,  particularly  elliptic  equations,  is  a  multilevel  technique,  to 
wit.  multigrid  (see  [Ref.  2],  [Ref.  3],  [Ref.  1]).  While  this  method  enjoys  remarkable 
success  in  solving  certain  classes  of  both  linear  and  non-linear  PDEs,  its  applicability 
to  solving  other  types  of  equations  awaits  development.  The  goal  of  this  work  is  two- 
fold: to  apply  the  FVE  method  to  discretize  a  particular  equation,  and  to  implement 
multigrid  in  solving  the  resulting  discretized  problem. 

The  approach  that  is  taken  is  to  begin  with  a  specific  physical  problem  which 
results  in  the  Navier-Stokes  equation,  and  consider  a  subsidiary  problem,  namely  the 
heat  equation;  due  to  the  nature  of  the  problem,  cylindrical  coordinates  are  used. 
The  resulting  problem,  to  which  the  FVE  discretization  and  multilevel  solution  are 
applied,  is  the  axisymmetric  heat  equation,  involving  the  use  of  a  point  source.  More 


specifically,  finite  differences  are  used  to  discretize  the  temporal  portion  of  the  prob- 
lem, resulting  in  the  use  of  a  fully  implicit  backward  time-stepping  scheme,  and  the 
FVE  method  is  used  to  discretize  the  spatial  part.  The  application  of  the  FVE 
met, hod  to  a  problem  in  cylindrical  coordinates  is  new.  and  results  in  stencils  which 
are  analyzed  extensively.  Once  the  discretization  has  been  constructed,  a  number  of 
relaxation  schemes,  including  both  Jacobi  and  Gauss-Seidel,  are  considered;  a  thor- 
ough analysis  of  these  schemes  is  done,  using  both  the  spectral  radii  of  the  iteration 
matrices  and  local  mode  analysis.  The  goal  is  to  choose  one  of  the  relaxation  schemes 
fur  use  in  an  iterative  solution  of  the  problem;  Gauss-Seidel  is  the  method  of  choice. 
The  specifics  of  the  multilevel  technique  are  then  developed  and  analyzed.  Local 
mode  analysis  is  performed  on  the  components  of  the  amplification  matrix,  resulting 
in  the  two-level  convergence  factors  for  various  combinations  of  the  operators.  The 
multilevel  solution  process  is  implemented  by  using  multigrid  V-cycles;  the  iterative 
and  multilevel  results  are  compared  and  discussed  in  detail.  The  computational  sav- 
ings resulting  from  the  multilevel  process  are  then  discussed.  Insofar  as  multigrid  is 
quite  successful  in  a  number  of  different  contexts,  the  results  of  our  work  are  some- 
what disappointing  in  that  we  are  not  yet  able  to  achieve  the  same  level  of  success. 
On  a  more  positive  note,  this  work  applies  the  FVE  method  to  a  problem  in  cylindri- 
ral  coordinates,  and  makes  use  of  Maple  software  to  compute  the  necessary  integrals 
to  construct  the  discretization  (see  Chapter  III  and  Appendix  B).  Additionally,  the 
techniques  that  are  applied  do  result  in  a  solution  to  the  problem. 


II.         PROBLEM  STATEMENT 

A.      BACKGROUND 

Our  ultimate  goal  is  to  be  able  to  numerically  solve  the  Navier-Stokes  equation. 

_^  +  £.Vu=  — Vp  +  i/V2u,  (III ) 

ot  p 

where  u  is  the  velocity  vector,  t  is  the  time,  p  is  the  density,  p  is  the  pressure,  and  v 
is  the  kinematic  viscosity.  One  approach  to  such  solutions  is  to  consider  this  problem 
,is  a  collection  ol  smaller  problems.  By  solving  each  of  the  smaller  problems,  we 
ran  assemble  the  collection  of  solution  processes  to  solve  the  original  problem.  The 
purpose  of  this  work  is  to  locus  on  an  initial  step  in  the  eventual  construction  of  a 
numerical  solution  to  Equation  ILL 

One  way  to  subdivide  a  problem  such  as  this  into  smaller  problems  is  to 
consider  a  specific  physical  problem  that  gives  rise  to  the  Navier-Stokes  equation,  and 
then  isolate  the  different  aspects  of  that  problem.  To  that  end,  we  consider  a  semi- 
infinite  block  of  metal,  with  a  heat  source  applied  to  the  horizontal  surface;  above  the 
horizontal  surface  is  an  (inviscid)  nonconducting  gas.  The  result  is  a  pool  of  molten 
metal  with  a  free  surface,  surrounded  by  the  solid  portion  of  the  unmelted  block  (see 
Figure  1).  The  total  heat  flux  Q  is  constant,  and  far  away  the  solid  approaches  the 
uniform  cold  temperature,  Tc.  The  resulting  thermal  and  flow  fields  are  assumed  to 
be  axisymmetric  and  steady,  and  are  governed  by  conservation  of  energy  in  the  solid 
and  by  conservation  of  energy,  momentum,  and  mass  in  the  pool.  We  therefore  have 
the  following  system  of  equations: 

solid         :     —  =  kV2T  (II.2) 

ot 

8T 
liquid        :     — +  u-VT  =  kV2T  (II.3) 

ot 

du       .    , 1 _  „0_  /IT1, 

— +  u-Vu= — Vp  +  t/V2u  (II.4) 

ot  p 

V-u  =  0  (II.5) 


here  T  is  the  temperature,  k  is  the  thermal  diffusivity,  and  u.  t,  p,  p,  and  v  are  as 


.1  hove 


Q 


Figure  1.  The  molten-pool  problem. 


The  portion  of  the  above  problem  which  is  our  focus  is  the  conduction  of  heat 
in  the  solid,  governed  by  the  heat  equation  (Equation  II. 2),  to  which  the  following 
boundary  conditions  apply: 


surface     z  —  0 

axis     /•  =  0 
far     away     r,  z  — >•  oo 


dT 
fc—  =  -q(r) 
oz 

dT      n 

7T  =  ° 
or 

T^TC, 


(II.6) 

(II.7) 
(II.8) 


where  k  is  the  thermal  conductivity,  and  q(r)  is  the  imposed  surface  heat  flux  (large 
at  /•  =  0,  falling  off  to  zero  at  some  small  value  of  r,  such  that  /0°°  q(r)27rrdr  =  Q);  the 
condition  in  Equation  II. 7  results  from  the  axisymmetric  nature  of  the  problem.  Since 
the  ultimate  goal  is  to  solve  the  Navier-Stokes  equation  in  cylindrical  coordinates,  the 
heat  equation  will  also  be  solved  in  this  coordinate  system. (see  [Ref.  4]) 


B.       POINT-SOURCE  PROBLEM 

Our  goal  is  to  take  the  first  step  toward  solving  the  Navier-Stokes  equation 
I  Equation  II. 1)  that  arises  from  considering  the  molten-pool  problem  by  solving  the 
associated  heat  equation  (Equation  II. 2).  We  therefore  assume  cylindrical  geome 
i  iv  and  azimuthal  symmetry.  We  further  simplify  the  conditions  imposed  on  Equa- 
tion II. 2  by  considering  that  the  heat  source  applied  to  the  horizontal  surface  of  the 
block  is  a  [joint  source,  with  the  total  heat  flux  Q  =  1,  and  that  far  away  the  solid 
approaches  the  uniform  cold  temperature  of  Tc  =  0.  That  is,  the  imposed  heat  flux 
</|rl  =  (S{r)  is  the  Dirac  delta  function,  with  /0°°  q[r)27crdr  =  1. 

The  existence  of  the  solution  to  this  type  of  problem  is  well  established.  How- 
ever, as  is  the  case  in  general  when  Neumann  boundary  conditions  apply,  absent  a 
specified  initial  temperature  distribution,  the  solution  is  not  unique  (see  [Ref.  5]). 
While  we  are  interested  in  solving  the  point-source  problem  numerically  in  cylindrical 
coordinates,  it  has  spherical  symmetry,  so  it  can  be  solved  analytically  in  spherical 
coordinates  more  easily.  Therefore,  for  the  analytic  solution,  we  rewrite  Equation  II. 2 
in  spherically  symmetric  coordinates: 

dr       i    d  (d2ot\ 

m=KWm[RdR)'  (IL9) 


where  R  is  in  spherical  coordinates,  R  =  \Jr2  +  z2.  Beginning  with  an  initial  tem- 
perature distribution  of  T  =  0  everywhere  and  Q  =  1,  the  solution  to  (II. 9)  with 
boundary  conditions 

surface     z  =  0    :     fc—  =  -q(r)  =  -S(r)  (11.10) 

az 

DT 
axis     r  =  0     :     tt"  =  0  (11.11) 

or 

far     away     r,  z  — >  oo     :     T  — >  0,  (11.12) 

is  found  to  be 

r^(»  =  ^berfc(A)'  (IU3) 


where  erfc(x)  is  the  complementary  error  function. 

2     yx    _2 

erfc(x)  =  I  —  ert(x).      erf(x)  =  — ■=   /     e   s  ds. 

yJTT  Jo 

As  t  — >  oo,  the  steady  solution  becomes 


T(R)  =  -±-  (11.14; 

27r/c/i 


( for  a  derivation  of  this  result,  see  [Ref.  5]).  In  (7-,  z)  coordinates,  where  R  =  Vr2  +  ~2. 
these  become 


1  (  Jx1  +  z2  \ 

T(r.zJ)     =     erfc     V      L-        ,      or  (11.15) 

r(r,,j     = i==.  (11.16) 

For  computational  purposes,  the  idealized  problem  of  a  semi-infinite  solid 
is  truncated  to  a  finite  domain  in  cylindrical  coordinates  with  azimuthal  symme- 
try. At  the  far  boundaries  on  this  finite  domain,  homogeneous  Dirichlet  conditions 
are  imposed  to  approximate  the  conditions  in  the  unbounded  problem.  To  further 
simplify  computational  work,  we  assume  that  the  (r, z)  domain  is  the  unit  grid 
VI   =   [0,  l]x[0,  1]   €   7£'2;we  will  eventually  assume  k  =   1.     Thus,  the  equation  we 

wish  to  solve  is 

3T 

-x-  =  kV2T,  (11.17) 

at 

with  boundary  conditions 

surface     z  =  0     :     k^-  =  -Sir)  (11.18) 

oz 

axis     r  =  0     :     %-  =  0  (11.19) 

far  boundaries  r,z  =  l     :    T  =  0.  (11.20) 


III.  DISCRETIZATION 

In  order  to  solve  the  conduction  problem  numerically  (Equation  1 1.2.  with 
boundary  conditions  Equations  II. IS,  11.19,  and  11.20),  the  equation  representing  the 
continuum  problem  must,  be  discretized.  That  is.  the  values  of  the  unknown  function 
f  are  to  be  determined  from  a  set  of  equations  which  in  some  .sense  approximate 
Equation  II. 2.  Ultimately,  a  discrete  approximation  to  the  heat  equation  should  meet 
the  following  criteria: 

1  )    Provide  for  a  unique  solution  to  the  problem. 

'_')   The  solution  should  be  "close"  to  the  exact  solution  for  "sufficiently  small" 
grid  spacings. 

■T)   The  solution  should  be  "effectively  computable.'1 

The  significance  of  property  1)  is  obvious;  property  2)  relates  to  the  questions  of 
convergence  and  consistency  for  the  discretization  scheme.  Property  3)  relates  to: 
a)  the  amount  of  computational  work  required  to  solve  the  problem,  and  b)  the 
behavior  of  roundoff  errors  in  the  computed  solution.  The  growth  of  the  roundoff 
nror  is  related  to  the  notion  of  the  stability  of  the  discretization  scheme,  (see  [Ref. 

«1) 

Solution  of  the  heat  equation  (Equation  1 1.2)  requires  treatment  of  both  space 
and  time  derivatives.  The  Finite  Volume  Element  (FVE)  method  is  used  to  discretize 
the  spatial  portion  of  the  problem;  finite  differences  are  used  to  discretize  the  tem- 
poral portion.  The  approximation  properties  of  the  FVE  method  are  fundamentally 
different  from  those  associated  with  the  finite  difference  method.  Hence,  this  ap- 
proach treats  time  and  space  differently.  The  goal  in  applying  the  FVE  method  to 
produce  spatial  discretization  is  to  make  use  of  the  advantages  of  the  Finite  Volume 
(FV)  method,  its  ability  to  be  faithful  to  the  physics  in  general  and  conservation 
in  particular,  coupled  with  the  flexibility  of  the  finite  element  representation  of  the 
unknown  functions  (see  [Ref.     1]).    As  outlined  in  [Ref.     1],  the  basic  approach  is 


to  partition  the  spatial  grid  in  two  ways:  as  the  union  of  a  set  of  finite  elements, 
the  vertices  of  which  comprise  the  grid  points  on  which  the  unknowns  are  defined 
(see  Figure  2),  and  as  the  union  of  a  set  of  control  volumes  (see  Figure  3),  one  for 
each  grid  point  (the  boundary  points  require  separate  treatment,  as  will  be  discussed 
later).  The  elements  are  used  to  form  basis  functions,  a  linear  combination  of  which  is 
used  to  approximate  the  unknown  function.  Upon  substitution  of  the  approximation 
into  the  equation  to  be  solved,  integrals  over  each  control  volume  are  taken.  The 
integrals  enforce  conservation  on  each  control  volume,  and  therefore  the  partition 
enforces  conservation  on  the  entire  domain.  The  system  of  equations  generated  by 
integration  yields  the  discretization  of  the  problem. 
N 


N 


Figure  2.  Partitioning  of  problem  domain  cross-section  H  by  elements. 

The  process  of  discretizing  Equation  II. 2  in  both  space  and  time  is  complicated, 
the  application  of  the  FVE  method  to  the  space  derivatives  being  the  more  complex 
portion  of  this  process.     Therefore,  before  proceeding  to  a  discussion  of  the  FVE 
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Figure  3.  Partitioning  of  problem  domain  cross-section  fl  into  sub-volumes. 

method,  the  finite  difference  method  is  applied  to  the  time  derivative.  In  order  to 
facilitate  this  presentation,  we  present  a  one-dimensional  example  of  this  type  of 
discretization  (see  [Ref.  7]).  Consider  the  one-dimensional  version  of  Equation  1 1.2. 

dT     d2T 


dt  "  dx2 ' 
where  k  —  1.  One  finite- difference  approximation  to  Equation  III.  1  is 

g  "  I?  ' 

where  T  is  the  exact  solution  of  the  approximating  difference  equations, 

Xfc  =  kh  and  tn  =  gn,      (fc,  n  =  0, 1,2,  •  •  •). 


Ill- 1 ) 


Let  ?•  =  A,  and  this  can  be  written  as 


Tk,n+i  =  rTk-i,n  -  (1  -  2r)rfc,„  +  rTk+i,n. 


Thus,  we  can  compute  the  unknown  7\,n+1  at  the  (n  +  l)st  time  step  using  the 
known  values  from  the  nth  time  step,  giving  an  explicit  formula  for  determing  the 
unknowns. [Ref.  7] 

While  this  explicit  relation  provides  a  simple  means  to  compute  unknown 
values,  is  has  a  serious  drawback.  The  process  is  onlv  valid  for  0  <  -^  h,  or 
tj  <  4r,  and  therefore  the  time  step  St  =  <j  is  necessarily  small.  Moreover,  in  order 
to  achieve  reasonable  accuracy,  Ax  =  h  must  also  be  kept  small.  There  is,  however. 
,i  method  which  reduces  the  total  amount  of  calculation  and  is  valid  (i.e.,  convergent 
and  stable)  for  all  finite  values  of  r,  which  was  proposed  by  Crank  and  Nicholson 
(see  [Ref.  7]).  The  technique  is  to  consider  Equation  III.l  as  being  satisfied  at 
the  midpoint  {kh.{n  +  -)g}  and  replace  —7  by  the  average  of  its  finite-difference 
approximations  at  the  nth  and  (n  +  l)st  time  steps.  That  is,  the  equation 

fdT\  (d2T\ 


is  approximated  by 


7fc,7i+l    "-  7\,n  f  7fc— l,n+l   "  ~Tk,n+l   +  T'fc+l.n  +  l  7fc-l,n  _  2Tfc,n  +  Tk+l,n 


g  2  1  h2  h2 

which  gives 

-''^-i.n+i  +  (2  +  2r)Tfc,n+i  -  rTk+i,n+i  =  rTk-i,n  +  (2  -  2r)T/t,n  +  rTk+hn, 

where  r  =  -j^.  Therefore,  we  now  have  three  unknowns  at  the  (n  +  1  )st  time  step 
written  in  terms  of  three  known  values  at  the  nth  time  step.  Thus,  the  Crank- 
Nicholson  method  generates  a  set  of  equations  which  must  be  solved  simultaneously; 
it  is  an  implicit  method. 

Although  the  Crank-Nicholson  method  for  Equation  III.l  is  stable  for  all  pos- 
itive values  of  r  in  the  sense  that  all  errors  eventually  tend  to  zero  as  n  tends  to 
infinity,  large  values  of  r,  e.g.,  40,  can  introduce  oscillations  in  the  solution  (see  [Ref. 
7]).    Therefore,  a  more  general  finite-difference  approximation  to  Equation  III.l  is 
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found  by  using  a  weighted  average  of  the  finite-difference  approximations  of  -^r  at 

i  lie  n\\\  and  ( n  +  1  )st  time  steps,  which  is  eiven  by 

o     — +  ( 1  -a 


'/  V  'n1  J  V  /i-2 

where  I)  <  a  <  1  is  possible.  Note  that  a  =  0  gives  the  explicit  scheme,  a  =  j.  Crank- 
Nicholson,  and  a  =  1  a  fully  implicit  backward  time-difference  method.  The  scheme 
is  unconditionally  valid  (stable  and  convergent)  for  \  <  a  <  1,  but  for  0  <  a  <  ^  we 
must  have  ?■  <  1/(2  —  4a).  Thus,  for  example,  for  a  =  0,  r  =  -jfa  <  t,  or  c;  <  -y,  is  a 
requirement  to  guarantee  validity. [Hef.  7] 

We  now  apply  this  finite-difference  method  of  discretizing  the  time  derivative 
iu  Equation  1 1.2.  Following  the  work  done  above,  a  general  finite-difference  approxi- 
mation for  the  time  derivative  of 

C\rp 

-i-  =  kV2T  (111.2) 

at 

is  given  by 

— — ^  =  o.KV2Tn+1  +  (1  -  a)/cV2Tn,  (III.3) 

9 

where  g  is  the  time  step,  and  0  <  cv  <  1  allows  for  a  weighted  average  of  the  current 
and  future  time  steps.  The  current  time  step  is  designated  by  n,  n  +  1  is  the  next 
time  step;  the  value  of  a  determines  the  type  of  time  stepping.  However,  we  have 
determined  that  in  order  to  guarantee  the  validity  (stability  and  convergence)  of  the 
explicit  version  of  this  scheme  (a  —  0),  we  must  have  g  <  y.  We  consider  this 
value  of  g  to  be  too  small  to  be  of  computational  use  (which  is  verified  by  numerical 
experiments),  and  therefore  we  will  use  the  values  a  =  ^  and  a  =  1.  Thus,  we  have 
the  semi-discrete  relationship  between  the  current  and  subsequent  time  steps: 

(I  _  aKV2)Tn+l  =  (-  +  (1  -  a)/cV2)Tn.  (III.4) 

9  9 

This  relationship  is  used  to  establish  a  discrete  set  of  equations,  which  gives  rise  to 

the  discretization  of  the  problem. 
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A.      BASIS  FUNCTIONS:  FINITE  ELEMENTS 

We  begin  the  discussion  of  discretizing  Equation  II. 2  in  the  spatial  variables  by 
recognizing  that  it  is  necessary  to  make  an  assumption  about  what  types  of  functions 
may  be  used  to  approximate  the  unknown  function  T:  we  consider  that  T  can  be 
approximated  by  continuous,  piecewise  linear  functions.  A  basis  is  then  chosen  for 
the  space  in  which  these  functions  are  found.  That  is,  a  set  of  functions  is  selected 
whose  linear  combinations  determine  the  space  of  functions  of  interest.  Therefore,  we 
choose  basis  functions.  <pk,j{ri  z),  which  are  assumed  to  be  continuous  and  piecewise 
lineai.  with 

^-)-fi;te,(r:).  (III.5) 

fc=0  ]=0 

The  next  step  is  to  determine  how  to  construct  these  basis  functions;  the  element 
partition  of  the  domain  Q,  is  used  as  a  framework  for  this  purpose,  (see  Figure  2). 

Since  we  are  assuming  cylindrical  geometry  with  azimuthal  symmetry,  the 
directions  of  interest  are  7'  and  z:  a  uniform  grid  of  step  size  h  =  jj  (N  =  number 
of  grid  points)  is  applied  in  both  directions.  In  this  way,  we  can  designate  function 
values  on  the  grid  Q  by 

Tk,j  =  T(rj,Zk), 

where  zk  =  kh,  r,  =  jh,  and  k,j  are  the  row  and  column  indices.  Each  square 
is  subdivided  into  two  triangular  elements  by  a  diagonal  (oriented  in  the  direction 
of  increasing  r  +  z)1;  the  basis  functions  to  be  used  in  approximating  the  unknown 
function  are  constructed  using  these  triangular  elements.  For  each  interior  grid  point, 
or  node2,  there  are  six  associated  triangular  elements:  NNE,  ENE,  SE,  SSW,WSW, 
NW  (see  Figure  4). 


'The  orientation  of  the  elements  is  based  on  the  necessity  to  apply  this  technique  to  solving  the 
molten-pool  problem.  In  particular,  we  will  be  required  to  track  the  movement  of  the  phase-change 
boundary  between  solid  and  liquid;  element  orientation  has  been  chosen  to  facilitate  keeping  track 
of  this  moving  boundary. 

2For  a  two-dimensional  grid  with  N  intervals  in  each  dimension,  there  are  (TV  —  l)2  interior  grid 
points,  (N  +  l)2  total  nodes. 
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Figure  4.    Hat  function,  top  view.   A:  is  the  row  index  (z  direction),  j  is  the  column 
index  [r  direction),  <p  direction  is  into  the  plane  of  the  paper  (0). 


If  we  suppose  that  this  collection  of  elements  is  joined  along  shared  borders, 
and  that  the  elements  can  be  "stretched"',  then  the  finite  element  for  a  particular 
node,  say  {k,j),  is  formed  by  anchoring  the  collection  along  its  external  boundary  (in 
the  (r,  z)  plane)  and  extruding  the  grid  point  in  a  direction  perpendicular  to  the  (r,  z) 
plane  and  parallel  to  the  tangent  to  the  ip  direction  at  the  grid  point  (k,j),  forming 
the  familiar  "hat"  function  (see  Figure  5).   Note  that 

1     at  the  node  (rj,  Zk) 

0      at  all  other  nodes 
and  is  piecewise  linear  over  each  triangular  element,  giving  the  hat  shape. 

By  choosing  these  hat  functions  to  be  our  basis  functions,  the  coefficients  of 
the  dk,j{r,z)  in  Equation  III. 5  become  the  values  /?jtj  =  Tk,j.  That  is,  we  can  now 
specify  the  nodal  values  of  T  on  the  grid  as  Tk,j  =  T(rj,Zk),  with 


<t>k,j 


Jt=0  j=0 


(III.6) 
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Figure  5.  Hat  function,  oblique  view. 

Since  the  4>k,]  are  piecewise  linear,  a  continuous  representation  for  function  values  of 
T  is  found  by  linear  interpolation  of  nodal  values  between  nodes.  For  example, 

Tk,j  +  7jfc,j+i 


Tk,>+k  = 


2 


Up  to  now,  we  have  considered  the  sampled  values  of  the  unknown  T  at  the 
(iY  +  l)2  grid  points  as  an  (N  +  l)x(N  4-  1)  two-dimensional  array;  it  will  later 
be  convenient  to  consider  T  as  a  one  dimensional  array.  One  simple  method  to 
accomplish  this  is  to  order  the  elements  Tk,j  lexicographically.  For  example, 

Tq,q  7o,i        •  •  •       TotN+\ 


\  Tn+ifl    Tyv+i,i     •••    7jv+i,ah-i  / 


can  be  written  as 


I  rp  rp  rp  rp  rp  rp  rp  rp  rp  \T 

U0,0,  io,l,  •  •  •  ,  io,N+l,  -11,0,  J-  1,1,  •  •  •  ,  il,AT+l,  •  •  •  ,  JjV+1,0,  ^JV+1,1,  •  •  •  ,  J-N+l,N+l)     , 
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which  we  will  normally  consider  to  be  a  column  vector.    Thus,  we  can  relabel  the 
.ihics  Tk.j  =  Ti,  I  =  ( N  -f  \)k  +  j  +  1  and  k.j  =  0  :  N  +  1,  in  Equation  III. 6  so  that 

T(r.  =  )=     £     T&iir,:).  (IIL7) 

B.       FVE  STENCILS 

With  the  construction  of  the  basis  functions  over  the  finite  elements,  we  can 
now  incorporate  the  control  volumes  to  complete  the  finite  volume  element  discretiza- 
t  ion.  Below  is  an  example  the  FVE  method  applied  to  a  simpler  problem:  the  purpose 
is  to  motivate  and  explain  the  work  that  is  necessary  to  handle  the  specifics  ot  apply- 
ing the  FVE  method  to  Equation  II. 2.  This  explanation  closely  follows  the  example 
in  [Ref.  1]. 

1.        Example:  2  — D  Potential  Flow 

The  example  problem  we  consider  is  the  potential  flow  problem  in  Euclidean 
two-dimensional  coordinates.  The  domain  for  the  problem  is  the  unit  grid  fi  = 
[0.  l]x[0,  1]  €  7£2,  where  x  =  0  and  y  =  0  are  the  near  boundaries  and  x  =  1  and 
y  =  1  are  the  far  boundaries.  The  governing  equation  is 

V.(/>Vtf)  =  i7,  (III.8) 

where  p  is  the  density,  xp  is  the  flow  potential  (i.e..  ipT  and  tpy  are  the  components 
of  the  fluid  velocity  in  the  x  and  y  directions),  and  7/  is  the  interior  source  flow  rate. 
The  boundary  conditions  are 

near     boundaries     :     (pVip)  ■  h  =  t/>i  (Neumann  boundary  condition), 
far     boundaries     :    x\)  =  ip0  (Dirichlet  boundary  condition), 

where  h  is  the  outward  normal.  Suppose  that  V  is  one  of  the  control  volumes  in  J7, 
similar  to  those  depicted  in  Figure  3,  where  r  and  z  are  replaced  by  x  and  y.  Since 
this  problem  is  in  two  dimensions,  the  control  volumes  V  are  actually  areas,  and  the 
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corresponding  surfaces  S  are  the  perimeters  of  these  areas.   Integrating  Equation  III. 8 
over  I   .  we  have 

/  V-(pVijj)dV  =   I  ridV. 

Jv  Jv 

By  applying  the  Gauss   Divergence  Theorem,  the  volume  integral  on  the  left-hand 
side  becomes  a  surface  integral. 

f  (pVtp)  •  hdS  =   fijdV.  (III.9) 

In  terms  of  the  physical  units  associated  with  the  potential  flow,  the  relationship  is 

mass       length  mass 

area  = ; •  volume; 


volume      time  time -volume 

each  side  represents  a  flow  rate  of  mass  per  unit  time,  and  (pVib)  ■  h  represents  a  flux 
across  S.  Therefore,  the  net  flow  from  the  interior  source  in  V  is  balanced  by  the 
flow  across  the  surface  S,  which  can  be  interpreted  as  a  mass  conservation  law  for  the 
control  volume  V'.[Ref.   1] 

By  partitioning  the  domain  ft  into  a  finite  collection  of  control  volumes  (see 
Figure  3,  where  r  and  z  are  replaced  by  x  and  y),  we  can  impose  on  each  control 
volume  the  condition  represented  in  Equation  III. 9  (the  boundaries  require  special 
treatment;  a  discussion  of  these  considerations  is  included  later).  As  noted  earlier, 
the  imposition  of  conservation  on  each  control  volume  results  in  conservation  over 
the  entire  domain.  This  integral  condition  on  each  control  volume  will  be  used  to 
compute  the  discretization;  the  number  of  discrete  equations  is  the  number  of  control 
volumes,  say  m,  used  to  partition  ft.  To  complete  the  discretization  process,  the 
unknown  function  ip  will  be  approximated  by  (3  in  7£m .  By  appropriately  choosing 
basis  functions,  we  can  construct  an  approximation  v  expressed  in  terms  of  its  nodal 
values.  Consider  a  triangular  element  partition  similar  to  that  in  Figure  2  (where  7' 
and  z  are  replaced  by  x  and  y)  and  let  T  be  the  space  of  continuous,  piecewise  linear 
functions  on  ft  associated  with  this  partition.  Ignoring  for  the  moment  conditions  at 
the  boundaries,  the  FVE  discretization  of  Equation  III. 8  comes  from  finding  auGT 
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in  ii  that  Equation  III.!)  holds  for  each  of  the  control  volumes  in  the  partition  of  il. 
lie  problem  in  7Zm  is  then  defined  when  v  is  expressed  in  terms  of  nodal  basis  for    f: 


m 


y=  E  "«*«»  (111.10) 

iu=i 

where  //„,  is  the  value  of  v  at  the  iyth  node  and  dw  is  the  "hat  function"  associated  with 

the  i«th  node  (as  above,  we  lexicographically  order  the  nodes,  resulting  in  a  single 

subscript).     Substituting   Equation   III. 10  into  Equation   III. 9,   we  have  the  matrix 

equal  ion 

Lu  =  /,  (III. 11) 

where  L  is  the  nxn  matrix  with  entries 

L^z  =  I    {pVo:)-MS  (III. 12) 

and 

U=  f    rjdV  (111.13) 

(except  for  boundary  conditions). [Ref.  1] 

The  treatment  of  the  boundaries  depends  the  type  of  boundary  condition. 
Neumann  conditions  can  be  imposed  indirectly  by  substituting  the  flux  value  ipi  into 
the  appropriate  term  in  Equation  III. 9.  Specifically,  suppose  we  choose  the  control 
volume  V  that  includes  the  origin  (lower  left-hand  corner,  Figure  3),  where  the  surface 
of  the  control  volume  coincides  with  the  boundary  along  the  two  axes.  The  integral 
condition  for  V  is 

/        i))XdS  +  /       (pVxb)  ■  hdS  =  [  rjdV;  (111.14) 

Jx,y=0  Jx,y^0  JV 

Equation  III.  14  is  substituted  for  the  interior  condition  in  Equation  III. 9  into  the 
discrete  approximation  v  for  the  corner  volume  containing  the  origin. [Ref.  1] 

Dirichlet  conditions  are  imposed  directly;  uw  takes  on  the  value  of  ip\  at  each 
Dirichlet  node,  i.e.,  each  node  on  the  far  boundaries  (x,y  =  1).  This  approach  may 
lead  to  fewer  unknowns  uw  than  equations,  a  problem  easily  resolved  by  discarding  the 
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equations  associated  with  the  Dirichlet  nodes.  In  this  case,  the  collection  of  control 
volumes  no  longer  partitions  the  domain  and.  therefore,  conservation  no  longer  is 
strictly  enforced.  However,  it  is  at  the  Dirichlet  nodes  where  the  loss  of  conservation 
is  not  a  concern  in  general. [Ref.   1] 

Assuming  that  there  is  a  balance  in  the  number  of  equations  and  unknowns,  we 
now  must  implement  a  numerical  rule  for  evaluating  the  integrals  in  Equations  III.  12 
and  III. 13.  For  Equation  III. 12.  we  use  the  following  rule  on  each  segment  A  that  is 
part  of  the  interior  surface  .5'  associated  with  a  volume  V: 

[  {pVS)  ■  hdS  ~  p{  PA )  [  Vo  ■  hdA 

■  ' A  J  A 

where  P\  is  point  of  intersection  of  A  and  the  grid  lines  passing  through  the  node 
associated  with  V .  For  Neumann  boundary  segments,  we  use  the  quadrature  rule 

/  iihdS~v,(PA)\Al 

J  A 

where  |.4|  is  the  "surface  area",  i.e.,  length,  of  A.  For  Equation  III. 13,  we  use 

/  ridV~7i(Nv)\V\, 

Jv 

where  Ny  is  the  node  associated  with  V  and  |V|  =  fydV  is  the  "volume",  i.e.,  area, 
of  \/.[Ref.  1] 

Numerically  evaluating  Equations  III. 12  and  III. 13  yields  the  FVE  discretiza- 
tion of  Equation  III. 8.  The  value  associated  with  each  node,  say  (p,  q),  is  a  sum 

w=p+l,z=q  +  \ 

/     „  ^UI,Z^W,Z  1 

w=p—  1  ,z=q—  1 

where  the  coefficients  aWtZ  are  determined  by  the  integration.  A  convenient  way  to 
represent  the  sum  associated  with  each  grid  point  is  to  place  the  coefficients  aWiZ  into 
a  stencil, 


aP+i,g 

&p,q  —  1  ®P,q  ^p,q+l 

OCp-l,q 


Uptq       —       OCp—\^qlLp^\q  -\-  Otptq  —  \1Lpyq—\ 

+       <Xp,qUp,q  +  OCp+ltqUp+itq  +  aP)(7+i  UP)(j+l 
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where,  for  this  example,  the  corner  values  in  the  stencil  are  all  zero.  The  value  for 
node  ip.ii)  ls  determined  by  applying  the  stencil  to  u.,  r 

To  see  what  type  of  stencils  are  produced,  consider  the  discretization  on  a 
uniform  mesh  with  grid  size  h  =  —  in  both  coordinates.    We  use  double  subscripts 

~  in  > 

]>.<l  =  I)  :  m  —  1.  where  0  corresponds  to  the  Neumann  boundary  nodes  and  in  —  1 
corresponds  to  the  Dirichlet  boundary  nodes.  Written  in  stencil  form,  the  interior 
nodes.  0  <  p,q  <  in  —  1.  for  the  equations  in  Equation  III. 11  are  as  follows  (J2 
indicates  the  sum  of  the  outer  stencil  entries): 


p{{p  +  \)k,qli) 
p[ph.{q-\)h)  -E  p{ph,(q  +  \)h) 

p((p-±)h,qh) 
For  the  corner  Neumann  node  (p,q  =  0),  the  stencil  is  given  by 


uVn  =  h2r](ph.qh). 


Mf,o; 


-E     M°4) 


"o,o  =  jr/(0,0)  -  £  f  0!(O,  £)  +  0i(pO) 


For  both  stencils,  the  coefficient  of  the  unknown  to  which  the  stencil  is  applied  is 

-E. 

For  grid  points  adjacent  to  a  Dirichlet  node,  the  stencil  entry  reaching  to  the 
Dirichlet  node  value  is  moved  to  the  right-hand  side  as  a  coefficient  of  the  boundary 
data.  Otherwise,  the  equations  for  these  points  are  the  same  as  above.  For  example, 
at  the  point  (0  <  p  <  m  —  \,q  —  rn  —  l)we  have 

p((p  +  J)M-/0 

p(pk,(q-l)h)  -E 

p((p-i)h,l-k) 


Wp,m-i     =  h2T](ph,  1  -  h) 


-p(phA  -  -)t/>o(pM), 
where  £  is  the  sum  of  the  outer  stencil  entries  without  the  boundary  terms  removed, 

J2  =  P(ph,  1  -  2 )  +  P(Ph>  (9  "  5  W  +  ^  +  2)/l' *  ~  /l)  +  pUp  "  2)/l' l  "  /i}" 
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see  IRef.    1" 


2.        Application  to  Axisymmetric  Heat  Equation 

This  technique  is  now  used  to  compute  the  FVE  stencils  for  Equation  1 1.2.  by 
combining  the  finite  elements  (see  Figure  2)  with  the  control  volumes  (see  Figure  3). 
A  control  volume  for  the  /th  grid  point,  call  it  V/,  is  a  toroidal  prism  (see  Figure  6). 
It  is  generated  by  taking  the  two-dimensional  sub-volume  for  that  point  (in  the  (r.  z) 
plane),  and  sweeping  through  360°  in  the  ip  direction. 


r-»    y 


control 
volume 


Figure  6.  Toroidal  volume  for  the  conduction  problem. 

Integrating  Equation  III. 4  over  all  control  volumes  VJ,  where  V 
partitions  the  domain  fi,  we  have 


=  u 


("+D2  V[ 


J  {--  <*k.V2\  Tn+ldV  =  J  (-  +  (1  -  a)KV2)  TndV,  (111.15) 
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or,  upon  application  of  the  (lauss  Divergence  Theorem  the  volume  integral  of  the 
V2T  term  becomes  a  surface  integral,  and  we  have 

-  /   Tn+ldV  -  cm  I  VTn+1  •  hdS  =  -  /  TndV  +  (1  -  <x)k  I  VTn  •  fidS.    (III. 16) 
(]  Jv  JS  (]  Jv  Js 

where  h  is  the  outward  normal.    Substituting  the  expression  for  approximating  the 
unknown  function  T  (Equation  III. 7).  we  have 


E 


(J  Jv  JS 


hdS 


rpil+l 


(111.17) 


E 


-  /  #W  +  ( 1  -  a)«  /  V<i£  •  fidS 
a  Jv  Js 


Ll    1 


(J  JV  JS 

where  we  now  have  integrals  over  each  of  the  (N  +  I)2  control  volumes,  resulting  in 
a  set  of  (iV  +  l)2  equations.  This  set  of  equations  can  be  written  both  as  an  operator 
equation,  L[Tn+I]  =  /(Tn),  and  as  a  matrix  equation3,  M[fn+1]  =  /(fn),  where  the 
operator  L  and  the  matrix  M  are  given  by 


a/cV2 


M     =     -  [  4>?+ldV  -an  f  V<f>?+1  ■  hdS, 
g  Jv  Js 

f{Tn)  and  /(fn)  are  given  by 

f(Tn)     =     |f-  +  (l-a)/cV'2|rW, 

-  /  tfdV  +  (1  -  a)K  I  V<f>?  ■  MS 
q  Jv  Js 


(N+l)2 

f(Tn)    =     £ 
/=1 


T,n. 


Any  function  whose  coefficients  satisfy  the  resulting  set  of  discrete  equations  neces- 
sarily satisfies  the  conservation  law  over  any  volume  made  up  of  the  union  of  control 
volumes,  except  possibly  at  the  boundaries.  The  Neumann  conditions  at  r  =  0  and 
z  =  0  are  incorporated  indirectly  into  T,  by  substituting  the  (zero)  normal  derivative 


3The  operator  may  be  represented  as  a  continuous  operator,  L;  as  a  discrete  operator,  L   ;  or  as 
a  matrix,  M. 
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value  into  the  appropriate  term  in  Equation  III. 15.  The  Dirichlet  conditions  at  the 
l.n  boundaries  are  imposed  directly  on  T:  control  volumes  that  would  usually  be  as- 
sociated with  these  boundary  nodes  are  eliminated  (see  Figure  7).  Thus,  the  reduced 
collection  of  control  volumes  no  longer  partitions  the  grid  fi.  which  slightly  impairs 
conservation  in  the  discretization.  However,  since  we  require  the  temperature  to  fall 
oil  to  zero  at  these  points,  this  loss  of  conservation  is  not  a  concern. [Ref.   1] 

N 


H 1 1 i 1 

_, 1 1 1 

n j j , 

n 1 1 1 


N 


Figure  7.  The  problem  domain  cross-section  is  depicted  after  partitioning  into  sub- 
volumes.  Homogeneous  Dirichlet  boundary  conditions  obviate  the  necessity  to  define 
volumes  for  grid  points  at  the  far  boundaries. 

The  integrals  over  the  basis  functions  have  been  computed  using  Maple  soft- 
ware and  a  program  designed  by  David  Canright  (see  Appendix  B  and  [Ref.  8]).  In 
order  to  calculate  the  integrals,  we  must  be  able  to  represent  the  unknown  function 
over  the  basis  functions,  which  themselves  are  collections  of  triangular  planes.  There- 
fore, the  method  is  to  determine  first  the  function  for  the  plane  through  three  given 
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points  (the  cases  of  vertical  planes  and  three  points  collinear  are  not  consider*.!  i 
I  Ins  procedure  is  then  used  to  create  the  triangular  elements  which,  when  assem- 
bled, form  the  hat  function  (see  Figure  5).  For  a  given  grid  point,  say  (/>.'. j),  six  of 
these  triangular  planes  are  joined,  corresponding  to  the  six  triangles  surrounding  the 
point.  NNE,  ENE,  SE,  SSW.WSW.  NW  (see  Figure  4).  The  unknown  function  is 
then  interpolated  over  the  six  elements. 

Once  the  unknown  function  is  represented  by  the  combination  of  these  six 
interpolations,  we  can  use  Maple  to  compute  the  volume  integrals,  and  the  surface 
integrals  of  the  gradients,  in  Equation  III.  17.  The  results  of  these  integrals  provide 
i  he  coefficients  which  comprise  the  FVE  stencils. 

Thus,  we  have  the  following  stencils  for  the  volume  and  surface  integrals  re- 
spectively for  interior  grid  points: 


/,(£XXa,,)^  =  |1 


'*■■.;       p       q 


32  j  -  5     16j  +  5 
32 j  -  1 1       224 j      32j  +  1 1 
16j  -  5     32 j  +  5 


TkJl         (111.18) 


and 


js  (£Evtp,^ms  =  ^ 


'fc,;       p       q 


V 


Tk,j  > 


(111.19) 


2j-l     -Sj    2j  +  l 

2j 

where  the  27T  resulting  from  integration  in  the  tp  direction  has  been  factored  out. 

(Recall  that  in  cylindrical  coordinates, 

I  dV  =  \U  f M   f k+*  rdzdrdtp.) 

JV  J0        Jr       ,       Jz.      , 

Note  that  control  volumes  increase  with  radial  distance  from  the  origin,  which  is 
reflected  in  a  radial  bias  in  the  stencils.  More  specifically,  Tj  =  jh  is  the  radial 
distance  from  the  origin,  with  j  the  column  index  for  the  unknown  matrix,  h  the 
meshsize.  Thus,  stencil  values  increase  with  distance  from  the  origin. 

These  stencils  are  applied  to  the  unknown  at  a  point  on  the  grid,  designated 
by  its  row/column  position  (k,j).   Stencil  entries  indicate  the  values  to  be  applied; 
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blank  entries  indicate  a  value  of  zero.  Entry  position  indicates  to  which  grid  point 
i  lie  value  is  applied:  left/right  and  up/down  in  the  stencil  correspond  to  neighboring 
points  in  those  directions  on  the  grid.  That  is.  if  stencil  values  are  replaced  by  the 
positions  to  which  they  apply,  we  have 

(kj-l)  (k,j)  (k.j  +  \) 

{k-l.j-l)     (k-l.j) 
rims,  lor  example,  the  value  that  appears  in  the  (k  +  1, j)  stencil  position  is  applied 
to  the  grid  point  that  occupies  that  same  position. 

Using  techniques  similar  to  those  outlined  in  the  two-dimensional  example, 
boundary  point  stencils  have  been  computed  for  the  volume  and  surface  integrals 
respectively.  At  the  origin,  7'  =  0,  z  —  0: 


^1 

384 

at  r  -  0: 


1      5 
15     3 


h 
Too,      and      - 

o 


1 

-3    2 


/«< 


To,o      +  /  q(r)rdr;         (111.20) 


384 


1      5 
25     11 
6 


Tk,o,      and 


1 
-6    4 

1 


Tk,o\ 


(111.21) 


and  at  z  =  0: 

32J-5     16J  +  5 
24.7-8     112j+5     8j  +  3 


384 


Tq  j  j      and      — 

4 


4j 
2j-l     -8j    2j  +  l 


^o,r 


(111.22) 

For  points  adjacent  to  the  far  boundaries,  the  stencils  in  Equations  III. 18,  III.  19,  III. 21, 
and  1 1 1.22  are  applied.  However,  since  the  homogeneous  Dirichlet  conditions  dictate 
that  these  boundary  values  are  zero,  the  resulting  contribution  after  stencil  appli- 
cation remains  zero.  Thus,  in  effect,  far-boundary  values  do  not  contribute  to  the 
stencils  for  points  adjacent  to  these  boundaries. 
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We  ran  now  combine  all  of  the  above  stencils  to  generate  the  operator  Lh , 


Lll[Tn+]]  =  /,l(Tn), 


11.23) 


where  h  is  the  step  size  on  the  grid.  The  matrix  representation  of  Lh.  M.  is  a  block 
fridia.gonal  (N  +  l)2x(N  +  l)2  matrix  of  the  form 


M  = 


A 

A 

0 

i) 

i) 

I) 

r 

.1 

A 

0 

0 

0 

i) 

r 

.1 

A 

0 

0 

0 

i) 

i) 

0 

i) 

0 

r 

A 

A 

0 

i) 

V 

A 

where  .4,  A,  and  T  are  ( N  +  1  )x(/V  +  1 )  generic  banded  matrices  (not  all  identical);  A 
is  tridiagonal.  A  is  upper  bidiagonal,  and  T  is  lower  bidiagonal.  When  M  multiplies 
the  matrix  of  (iV  +  l)2  values  of  T,  arranged  lexicographically  as  a  column  vector, 
the  matrices  A,  A,  and  V  produce  the  effect  of  the  stencils  "reaching"  function  val- 
ues respectively  on  the  current  row,  the  row  above,  and  the  row  below.  Numerical 
experiments  using  Afatlab  to  construct  the  matrix  M  for  iV  =  8,  12,  16,20,  and  24, 
with  g  =  h  and  a  =  ^  and  1,  indicate  that  it  has  full  rank.  Thus,  we  expect  a  unique 
solution  to  the  linear  algebra  problem  which  is  a  discretization  of  the  heat  equation. 
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I IV 


IV.  RELAXATION  SCHEMES 

The   FVE  method,  with  weighted-average  time-stepping,  has  been   used   to 
iscretize  the  continuum  problem. 

—  =  kV27\  (iv. n 

<)t 

Lh[Tn+l]  =  fk{Tn)  (IV.2) 

or.  in  matrix  form. 

M[r,  +  I]  =/(fn),  (IV.3) 

on  a  grid  of  meshsize  h  =  y.  The  input  for  these  equations  is  the  solution  at  the 
current  time  step,  Tn  (where  Tn  and  Tn  are  used  interchangeably),  the  result  of 
solving  the  equations  is  the  solution  at  the  next  time  step,  Tn+1 .  As  was  indicated 
in  Chapter  III,  there  is  good  reason  to  believe  that  the  matrix  M  is  of  full  rank  and. 
thus,  expect  a  unique  solution  to  exist  for  the  linear  algebra  problem  that  arises  from 
the  discretization  of  the  continuum  problem.  Solution  by  direct  methods  requires 
factorization  of  M  and,  since  M  is  both  large  and  sparse,  solution  by  direct  methods 
may  1hj  impractical.  We  therefore  consider  iterative  methods  to  solve  the  matrix 
(■((nation  at  each  time  step.  These  methods  generate,  for  each  time  step,  a  sequence 
of  approximate  solutions,  (T(s) } l ;  the  choice  of  relaxation  scheme  determines  whether 
or  not  the  sequence  {T(s)}  converges.  Upon  implementation  of  a  relaxation  scheme 
that  produces  a  converging  sequence  of  approximate  solutions  at  each  time  step,  we 
begin  with  an  initial  temperature  distribution,  T°,  and  the  solutions  are  stepped  in 
time.  The  desired  result  is  a  sequence  of  solutions,  corresponding  to  a  sequence  of 
time  steps,  which  tends  to  steady  state.    This  process  allows  for  the  evaluation  of 


1  Here  the  subscript   (s)   indexes  the  sequence  of  approximate  solutions;  elsewhere,  subscripts 
without  parentheses  indicate  grid  position  or  vector  components. 
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various  values  for  the  time  step,  as  well  as  of  various  weightings  used  in  averaging  the 
current  and  subsequent  time  steps. 

A.      ITERATION  MATRICES 

It  is  common  in  constructing  iterative  methods  to  propose  a  splitting  for  the 
matrix  M  =  E  —  F.  where  linear  systems  of  the  form  Ex  =  b  are  ''easy"  to  solve  (see 
[R.e.f.  !)]).  The  sequence  of  approximations.  {T^},  is  generated  by  ET(j+l|  =  FT(sj+/ 
and.  therefore,  it  is  natural  to  define  an  iteration  matrix.  P  =  E-1F.  so  that 
I^s+i)  —  P^(s)  +  E-1/-  Additionally,  if  T*  is  the  exact  solution  so  that  WIT"  —  /, 
then  ETm  =  FT"  +  /  and  f*  =  E~lFT*  +  E"1/.  Hence  the  solution.  f".  is  a  fixed 
point  of  the  iteration  Tt3+n  =  E_1FT(5)  +  E-1/. 

The  matrix  P  =  E~'F  is  also  called  the  error  propagation  matrix,  since  if 
T(s+\)  =  P^(5)  +  E-1/  or  T(s+i)  =  PT(sj  +  B  (B  a  constant  vector),  and  e(s+1)  is  the 
error  at  the  (s  -f  1  )st  step,  then 

e(s+i)     =    ^(s+i)  -  T* 

=    [Pf(s)  +  B)  -  [Pf *  +  B]  (since  f*  =  Pf*  +  B) 
=    PTW-P!P 

=  p(f(s)-n, 

and  therefore 

e(s+i)  =  Pe(5).  (IV.4) 

Using  induction,  it  is  easy  to  show  that  e*(s)  =  Pse*(0).  By  Equation  IV.4,  e^  =  Pe(0) 
and 

e{2)    =    Pe(i) 

=    P[Pc(o)] 
=    P2e(0). 
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By  induction,  assume  that  ea.)  —  Pkff()),  in  order  to  show  that  e^+i)  =  P(^+"f(,M. 
Now. 


e(k+\ 


=     Pc(*, 

=       P[P%, 


j  (by  the  induction  hypothesis) 


_      p(*+Dr 


=   p 


r(o)- 


Thus,  since  G(k+i)  —  Pl*:+1)£r(o)i  we  conclude  e*(s)  =  Psc*(0),  for  s  any  integer.    Addi- 
i  ionallv. 


|es||  =  ||P*c(0,||  <  ||P 


'(0)1 


IV.51 


which  will  be  useful  later. 

One  of  the  simplest  relaxation  methods  is  the  Jacobi  method,  also  referred 
to  as  simultaneous  displacement  (for  a  more  detailed  account  of  all  of  these  methods, 
see  [Kef.  '.]}  or  [Ref.  9]).  It  is  produced  by  solving  the  /th  equation  of  Equation  IV. 3 
for  the  /th  unknown.  T),  /  =  1  :  (N  +  I)2.  Before  proceeding  to  a  discussion  of  the 
iteration  matrix,  we  present  the  component  form  for  this  iteration  scheme: 


fcj 


=  77 


/i-  [^(EvV^er 


_|i.224j  +  ^8j 


(iv.e: 


where  £^y  and  ]T5-  are-  respectively,  the  sum  of  volume  stencil  entries  and  surface 
stencil  entries  applied  to  the  current  values  of  the  unknown,  T£/d,  except  at  the  grid 
point  on  which  the  stencil  is  centered: 
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-    old 


Ev  = 


,(old) 


(old)  \ 


(32j  -  5)2^      +{lf>j  +  b)T£?J+l 


+(32j--ll)rg?, 

,(oM) 


(o/d) 


+(32J  +  11)T^;'1 


-(oW) 


and 


^ +(i6i  -  5)rrri_,  +(32j  +  5)rni 


(IV.7) 


-  yold 


E,  = 


(old) 


^■rj-^Oia, 


(oW) 


+(-1  +  2i)r£:'1 


+(i  +  2i)r« 


(oW) 


+yrni 


(IV.8) 
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One  iteration  sweep  consists  of  computing  Equation  IV. 6  for  each  component  of  the 
unknown  vector  T .  Provided  that  the  process  converges,  the  sweeps  continue  until  a 
desired  level  of  convergence  is  reached. 

Another,  more  succinct,  way  to  present  this  method  is  in  matrix  form.  If 
we  write  the  operator  matrix  as  the  sum  of  a  diagonal  matrix  (D),  a  strictly  upper 
triangular  matrix  (U),  and  a  stictly  lower  triangular  matrix  (L), 

M  =  D  +  U  +  L, 

the  matrix  equation  to  be  solved  becomes 

(D  +  U  +  L)[f]  =  /. 

Now  define  E,/  =  D  and  Fj  =  —  (U  +  L),  so  that  this  may  be  written  as 

D[f]     -     _(U  +  L)[f]  +  /,     or 
Ej[f]    =    Fj[f]  +  f, 

or  as 

f    =     -D-1(U  +  L)[f]  +  D"1/,     or 
f    =    E?Fj[f]  +  Ejlf, 

which  corresponds  to  solving  the  /th  equation  for  7/,  /  =  1  :  (N  +  L)  ■  If  we  define 
the  Jacobi  iteration  matrix  as 

Pj     =    Ej'Fj,      or 

=    -D-J(U  +  L), 

the  matrix  form  of  the  Jacobi  method  becomes 

f(new)   _  Pj[f  (oW)j  +  E-lf  (IV.9) 
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A  modified  version  of  this  method,  called  the  weighted  Jacobi  method,  is 
iletermined  by  introducing  a  weight.  U  <  uJ  <  1  (^j  =  1  is  the  original  .Jacobi  method) 
(see  [lief.  •!]).    The  matrix  form  is 


T(new)  _  p^T{oid)}  +cc,e-'/. 


IV.  10' 


where  I  is  the  indentity  matrix.  E^  =  D.  and  F^  =  (1  —  lj)D  —  u;(U  +  L).  and 

P      =    E~'F 

=    D"l[(l  -u;)D-lj(U  +  L)] 

=      (1  -i^I  +  ^Pj 

is  the  weighted  Jacobi  iteration  matrix.  In  this  method,  the  new  approximation  is  a 
weighted  average  of  an  intermediate  approximation  and  the  old  approximation;  the 
intermediate  approximation  is  the  Jacobi  iterate  of  the  old  approximation. 

The  weighted  Jacobi  method  computes  all  of  the  components  of  the  new  ap- 
proximation before  it  begins  to  use  them  in  the  next  approximation.  This  requires 
storing  both  the  current  and  new  approximations;  a  simple  modification  allows  for 
updating  the  current  approximation  "in  place",  using  only  the  storage  required  for 
one  approximation.  The  Gauss-Seidel  method  incorporates  the  new  changes  as  soon 
;is  they  are  computed  by  overwriting  the  current  approximation  with  the  new  (see 
[lief.  '.\\).  More  important,  however,  is  that  (on  model  problems)  the  Gauss-Seidel 
method  generally  converges  about  twice  as  fast  at  the  Jacobi  method  (see  [Ref.  9]). 
In  component  form,  this  iteration  scheme  is 


rp(new)   _  rp(new)  *l         1^384^^)  2     \£-S)\ 
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where  the  arrow  indicates  overwriting,  and 
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,111(1 
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(iv.i.r 
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\  +Vt£\1  ) 

11  matrix  form,  using  M  =  D  +  U  +  L.  Ery  =  (D  +  L),  and  Fq  =  — U,  we  have 

(D  +  L)[f]     =     -U[f]  +  /,      or 
EG[f]    =    FG[f]  +  /, 

f(new)     ^    _(D  +  L)-1U[f(o/d)]  +  (D  +  L)-1/,     or 

f(»««)    <_   e^FcCT*0'10]  +  E51/. 

This  corresponds  to  solving  the  /th  equation  for  7),  using  the  new  approximations  for 
components  1,2, ... ,/  —  1.  Now,  define  the  Gauss-Seidel  iteration  matrix, 

PG    =    EalFa 

=    -(D  +  L)-]U, 


so  that  the  iteration  scheme  in  matrix  form  becomes 


f(new)  j_  pG[j(oW)j  +  E"1/. 


(IV.14) 


With  the  above  lexicographic  ordering  of  the  (/V  +  l)2  components  of  T  and  the 
components  updated  in  ascending  order,  the  effect  of  a  Gauss-Seidel  sweep  is  to  start 
at  r  =  0  and  update  in  the  radial  direction  for  each  vertical  step,  starting  at  z  =  0 
(see  Figure  8). 

As  with  the  weighted  Jacobi  method,  we  can  make  a  modification  to  the  Gauss- 
Seidel  method.  Define  a  parameter  7  E  71  (7  =  1  is  the  original  Gauss-Seidel  method) 
and  the  modified  method  is  successive  over-relaxation,  SOR,  (see  [Ref.  9])  which, 
in  matrix  form,  is  given  by 


f{neu>)  =p^f(^)+7E"1/, 


(IV.15) 
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0  old  approximation 

1  new  approximation 
Figure  8.  Gauss-Seidel  sweep. 


where,  with  E7  =  D  +  7L  and  F7  =  (1  —  7)D  —  7U,  we  have 

P7    =    (D  +  7L)-,[(1-7)D-7U] 

=   e;'f7. 

Similar  to  the  weighted  Jacobi  method,  SOR  is  a  weighted  average  of  an  intermediate 
approximation  and  the  old  approximation. 

The  question  of  interest  now  is  whether  or  not  the  sequence  of  approximations. 
{T(4.)},  generated  by  ET(5+1)  =  FT^j  +  /,  converges.  Therefore,  we  make  use  of  the 
following  theorem: 

Theorem  IV.  1  Suppose  that  f  G  7ln  and  M  =  E  -  F  G  7lnxn  is  nonsingular.  If 
E  is  noTisingular  and  the  spectral  radius  of  E-1F  satisfies  p(K~lF)  <  1,  then  the 
iterates  T(s),  defined  by  ET(s+1j  =  FT(5)  +  /,  converge  to  T  =  M-1/  for  any  starting 
vector  T(0). 

The  proof  is  found  in  [Ref.  9],  and  makes  use  of  the  following  lemma: 

Lemma  IV. 1   //p(E_1F)  <  1,  then  (E_1F)*  ->■  0  as  k  ->  00. 
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N  =  8  |  .V  =  12 

N  =   16 

N  =  20 

.V  =  24 

Jacobi 

.9259 

.9544    .9675 

.9749 

.9796 

w=  L 

.9721 

.9826 

.9875 

.9903 

.9921 

"=  i 

.9442 

.9652 

.9750 

.9806 

.9841 

Gauss-Seidel 

.8395 

.8982 

.9263 

.9425 

.9530 

7^ 

.9670 

.9793 

.9851 

.9884 

.9905 

7  =  ^ 

.9189 

.9488 

.9630 

.9711 

.9764 

Table  I.  Spectral  Radii  for  Crank-Nicholson  Time  Stepping  (a  =  j 


N  =  8 

X  =  12 

.V  =  16 

X  =  20 

X  =  24 

Jacobi 

.9499 

.9709 

.9801 

.9850 

.9881 

u,  =  l 

.9817 

.9892 

.9925 

.9943 

.9955 

s 

.9634 

.9784 

.9850 

.9887 

.9909 

Gauss-Seidel 

.8931 

.9362 

.9556 

.9663 

.9730 

7=^ 

.9782 

.9871 

.9910 

.9932 

.9946 

) 

7  =  5 

.9462 

.9680 

.9777 

.9831 

.9865 

Table  II.  Spectral  Radii  for  Implicit  Time  Stepping  (a  =  1). 

The  matrices  M,  Ej,  and  Eg,  and  iteration  matrices,  P,  for  various  values 
of  Ar,  a  =  ~2  and  1,  and  g  =  /i,  have  been  constructed  using  Matlab.  We  have 
experimentally  verified  that  M,  Ej,  and  Eg  are  nonsingular,  and  the  spectral  radii  of 
the  error  propagation  matrices,  P.  have  been  computed.  Due  to  memory  limitations 
with  Matlab,  the  calculations  are  made  for  relatively  small  values  of  N.  The  results 
are  indicated  in  Tables  I  and  II,  where  the  value  of  a  indicates  the  type  of  time 
stepping,  and  the  values  of  u  and  7  are  the  weightings  in  the  weighted  .Jacobi  and 
SOR  methods  respectively. 

The  results  of  these  numerical  experiments  indicate  that  both  p(Pj)  <  1  and 
p{Pa)  <  1,  with  p(Pg)  <  p(Pj).  The  modifications  to  the  Jacobi  and  Gauss-Seidel 
methods  do  not  provide  any  improvement;  the  spectral  radii  for  both  methods  are 
higher  for  the  modified  iteration  matrices  than  for  the  original  iteration  matrices. 
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Iliiis.  the  (Jauss-Seidel  method,  used  with  either  the  Crank-Nicholson  or  implicit 
I  in le-slepping  scheme,  appears  to  i>e  the  besl  choice  of  these  relaxation  schemes.  How- 
ever, il  the  spectral  radius  of  P  is  near  unity,  convergence  may  by  unacceptably  slow 

i  uce  the  error  tends  to  0  likep(P)  .  .is  indicated  by  the  lemma  and  \\rs  _  ||P||"||t(o)|| 
(Equation  IV. o).  As  is  evident  from  the  above  tables,  even  for  a  moderately  spaced 
grid  (e.g.,  .V  =  16,  </  =  h,  and  a  —  \  or  1),  p(P)  >  .9,  and  the  spectral  radius 
increases  with  S . 


B.       ALTERNATIVE  SCHEMES:  LINE  RELAXATION 

Both  the  Jacobi  and  Gauss-Seidel  methods  give  rise  to  iteration  matrices,  are 
implemented  in  a  straightforward  manner,  and  are  attractive  because  of  their  sim- 
plicity. While  the  ease  of  implementation  is  a  significant  advantage,  the  convergence 
properties  of  either  or  both  may  not  be  acceptable  and,  therefore,  alternative  schemes 
may  be  desirable.  One  such,  which  seems  reasonable  based  on  the  geometry  of  the 
problem,  is  line  relaxation.  There  are  two  obvious  options  in  this  regard:  radial  or 
vertical  line  relaxation  (see  Figures  9  and  10),  where  either  an  entire  row  or  column 
is  updated  simultaneously.  Both  options  require  the  solution  of  an  (N  +  l)x(;V  -f  1) 
tridiagonal  matrix  for  each  row/column  update  in  the  unknown  matrix.  That  is.  while 
the  previously  outlined  relaxation  schemes  (point  relaxation)  proceed  by  successively 
solving  a  collection  of  algebraic  equations  for  one  unknown,  line  relaxation  proceeds 
by  solving  a  succession  of  matrix  equations.  For  example,  suppose  that  the  matrices 
.4.  A.  and  T  are  tridiagonal,  upper  bidiagonal.  and  lower  bidiagonal,  respectively,  and 
that  T,  and  /,  are  the  portions  of  the  respective  unknown  and  right  hand  side  vectors 
in  Equation  IV. 3  that  are  rows  in  their  corresponding  matrices: 

A    A     0     0 
T     A    A     0 

o    r    .4    A 

o    o    r   a 


"r, " 

'  h 

T2 

h 

T3 

h 

.T«  . 

./«. 
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Radial  line  relaxation  consists  of  solving  the  following  succession  of  matrix  equations: 


T, 


A-H/2  -rr,-Ar3), 

A-l(f3-rT2-AT4),  and 


t4  <-  A-^A-rra). 


In  order  to  illustrate  how  line  relaxation  works  relative  to  the  FVE  stencils. 


let 


/ 


r  >iu 


(32j  -  5)7&'*       +(16j  +  5)j£ft  +1  ' 


i(7ieto) 


i(?ieti;) 


and 


^+(16^-5)7^,     +(32i  +  5)rrr3 

/ 

rotw 


v 


y-rp(old) 

+2jt£z] 


and 


column 

£  = 

v 


(neui) 


+(32j  -  ll)Tfc7_l 


(16i  +  5)Ti^.+1    ^ 

+(32j  +  n)r« 


netf) 


lv+(16J-5)Tri7-i 


and 


Co/um?l 


£  = 


/ 


v 


\ 


(new) 


+(-l+2j)T^ 


(o«) 


+(i  +  2j)r^ 


Now,  define 


/i3 


IV.16) 


(IV.17) 


(IV.18) 


(IV.19) 


(new) 


^(new)■\ 


[(32 j  -  U)rin  +  224jTi7""  +  (32  j  +  U)T^] 


S3841 
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( new)column 


1  k.j 


h[i 


-[(32j  +  o)T£Z]  +  224i^i7l',  +  (32j  -  o)T^'] 


</384 


so  that,  for  the  radial  and  vertical  relaxations  respectively,  we  must  solve 
radial     :     T<7"r<""  -  /,  -  [—-(£.)  "  ^(Dl     '" 


vertical 


y*(tiew)coiinnn 

1  fc,i 


V/384^ 

;  3  '  Wiurni 


a/c 


/j 


column 


J'-'W4{    £    ,"~(    ¥    " 


simultaneouslv  for  each  row  or  column  in  the  unknown  matrix. 


t 

t 

^ next  row  update 

^ latest  row  update 


(IV.20] 

(iv.2i; 


—  -  -       old  approximauon 
^^^™       new  approximation 

Figure  9.   Updating  an  entire  row  at  one  time. 

The  computational  cost  of  updating  an  entire  row/column  at  a  time  is  higher 
than  a  row/column  update  using  either  the  Jacobi  or  Gauss-Seidel  methods.  However, 
this  type  of  relaxation  may  be  sufficiently  efficacious  to  warrant  the  extra  expense,  in 
that  convergence  may  be  achieved  significantly  faster.  For  comparison,  we  now  con- 
sider local  mode  analysis  of  the  Jacobi  and  Gauss-Seidel  methods  and  line  relaxation. 


:*7 


latest  column  update  next  column  update 


—  ~  —        old  approximation 
^^^™        new  approx imation 

Figure  10.  Updating  an  entire  column  at  a  time. 

C.      LOCAL  MODE  ANALYSIS 

While  an  examination  of  the  spectral  radii  of  iteration  (error  progagation) 
matrices  can  be  instructive,  it  is  often  useful  to  conduct  a  more  detailed  analysis. 
\\  hat  follows  is  a  local  mode  analysis  of  various  iteration  schemes.  Since  the  error 
propagation  matrices  indicate  how  the  error  evolves  during  the  iteration  process,  we 
use  a  DFT  (see  Appendix  A,  [Ref.  10])  to  expand  components  of  the  error  equations 
associated  with  the  various  relaxation  schemes.  The  coefficients  in  these  expansions 
are  the  factors  by  which  the  corresponding  modes  of  the  error  are  magnified/reduced 
for  each  relaxation  sweep.  Thus,  by  examining  these  coefficients,  we  can  determine 
how  quickly  the  error  tends  to  zero,  which  indicates  how  quickly  the  solution  con- 
verges. The  following  analysis  does  not  apply  to  points  on  the  boundary,  it  is  only 
valid  for  interior  points.  Thus,  while  it  gives  more  information  about  the  functioning 
of  the  scheme  in  the  interior  of  the  domain,  boundary  peculiarities  are  not  addressed. 
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We  begin  by  recalling  that,  the  current  error.  e(s),  is  the  difference  between 
i  he  exact  solution.  T*.  and  the  current  approximation.  Tis)i 

e*(s)  =  T*  —  /(s). 

where  we  desire  that  the  sequence.  {f(s)}.  tend  to  zero.  Substituting  this  expression 
into  Equations  IV. (>  and  IV. 11,  we  have  the  following  error  equations  for  the  Jacobi 
.uid  Gauss-Seidel  methods  where,  in  order  to  avoid  confusion  with  the  exponential 
[unction  in  the  DFT  expansion,  we  represent  the  components  of  e  by  v^j- 


,  ,  •  (new)  j384V^V    I  j     VZ-.s     . 

lacobi     :     o\  ■       =  —  - 


-i£- 224  7  +  —  8/ 

h3    /v-      \         a/s/i. 


IV. 22' 


Ciauss  —  Seidel  yv  •       =  —  - — n t ,  IV. zJ 

k,j  _fcL_224?  +  — 87 

g384 *•>    ^      2     U^ 

where  Ylvi  Ylsi  Hvi  and  zZs  are  now  applied  to  the  Vk,j  (see  Equations  IV. 7,  IV. 8,  IV.  12, 
and  IV.  13).  Equations  IV. 22  and  IV. 23  indicate  how  the  sequence  {e*(s)}  evolves 
during  the  iteration  process.  If  we  expand  these  relationships  in  a  discrete  Fourier 
transform,  (DFT)  (see  Appendix  A  or  [Ref.  10]),  the  magnitude  of  the  transform  co- 
efficients will  indicate  the  "amount"  of  each  mode  of  the  error  that  is  present  in  each 
component  of  eu\.  We  then  compare  coefficients  to  determine  the  ratio  between  the 
new  and  current  approximations  for  each  component  of  the  error.  In  other  words,  the 
DFT  allows  us  to  analyze  the  growth  of  the  error  by  examining  component  behavior. 
Expanding  Equation  IV. 22  in  a  DFT,  where 


7-nkl       '2irjm 


\7nkl       i2ir;m 


1 


C(l,m)  =  e   n   e    n     -^  0,      and     p(j)  =  -73— 

and  the  superscripts  (0)  and  (n)  indicate  the  components  of  the  old  and  new  approx- 
imations, we  have 
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V  V 


'/.?n  =  — £  + 


V^CiLm)    = 


i 


h 


5-384 

2n(l+m) 


((32i-5)ViS;C(/,m)c^ 


+  (16j  +  5)V/£)C(f,m)e 
+(32j-ll)V^C(/,m)e^ 
+(32j  +  ll)V&)C(/,m)e 


-i2_irm 
.V 


t27rm 
N 


i2ir((+ml 


+(16j-5)^C(/,m)e 

+(32j+5)VSC(/,m)e-^) 
—  (2jV;(;i)C(/.m)eN 


+(-i  +  2j)v;;{;i)c(z,m)( 

+(l+2j)V^C(/,m)e^ 


:27rm 
.V 


-i2tt( 


+2jVSC(/,m)e^)]}J-(j), 


Making  use  of  the  orthogonality  property  of  the  complex  exponential  (see  Appendix  A) 
we  can  equate  individual  terms  of  the  sum,  and  then  divide  by  C(l,m)  to  give 


V, 


(n) 


l,m 


H  t„\      i2nl  (r,\      t2ir(l+m) 


,(fl 


(°)^ 


+(32j  -  llW^e-ir-  +  (32j  +  U)V^e 
+(16i  "  5)K£e-^^  +  (32,  +  5)l^e-V) 


anil 


'(°L%* 


(°)„-* 


—ri^Zle^r +(-l+2j)V^e 


+(i+»^+^'^)  na 


or 


((32j  -  5)e—  +  (16;  +  5)e"^^ 


2»rJ 


^384 
+(32j  -  l^e^21  +  (32 j  +  \\)e*% 

i2jr(l  +  m)  „  i„„ 

+(16.7  -  5)e TT^  +  (32j  +  5)e"— ) 

Qlft/i   .  i2-rrl  .  „    ..  i2ffm 

— (2je—  +  (-1  +  2j)c— 7T- 

+(1  +  2;^  +2jeldPj\ru)vQ 
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N=  16 

N  =  32 

i 

I   1    •       .v   ' 

=  ;V  -  !  ! 

j=  1 

_    A_ 

J  =  V  -  i  | 

1  tt=i 

L.009 

1.009 

1.009 

1.026 

1.026 

1.026 

1  a  =  1 

1.005 

1.005 

1.005 

■ 

1.013 

1.013 

1.013 

Table   III.  Maximum  Values  of 
Relaxation. 


r(n). 


'{«)\ 


for  /,  m   —  0   :    N  for  the  Jacobi   Method  of 


herefore.  the  ratio  of  the  new  coefficient  to  the  old  is 


An) 

vl.m 


k3  ,2rrl  .2»rH+r,.i 

((32j-5)e—  +(16j+5)c— ~ 


^384 


+(32j-  U)e^-  +  (32j  +  ll)e^~ 
+(16j  -  5)e V-^  +  (32j  +  5)e-~) 


a/c/ 


(2je  w    4-  (-1  +  2j)e      n 


+(1  +2j)e~^"  +2je 


Hj)- 


IV.24) 


In  older  to  determine  the  greatest  factor  by  which  a  mode  of  the  error  is  multiplied. 

we  seek  the  maximum  of     ';"%    over  the  values  /,  in  =  —  ^7  +  1  :  *%■.  In  other  words,  we 

ivgl  2  2 

seek  to  determine  a  bound  on  how  well  we  can  expect  this  type  of  relaxation  scheme 
to  perform.  This  ratio  is  a  function  of  NJ,tn  and  j  and,  while  difficult  to  determine 
analytically,  may  be  calculated  numerically.  Using  Matlab,  the  maximum  of  this  ratio 
for  /.  m  =  0  :  N  has  been  calculated  for  N  =  16,32;  a=U;  and  j  =  1,  y,  N—  1;  the 
results  are  indicated  in  Table  III  (see  Appendix  A  for  a  discussion  of  the  equivalence 
of  centered  indices,  l,m  =  — -y  +  1  :  y,  and  non-centered  indices,  /,m  =  0  :  :V). 
Additionally,  by  considering  the  matrix  of  grid  points  /,  m  =  0  :  N,  we  can  determine 

a  correspondence  between  sample  points  on  the  grid  and  type  of  associated  frequency 

|v(n)| 
(see  Figure  11  and  Appendix  A).  The  matrices  of  values  of     '(™    for  j  =  N  —  1  are 

l*J,ml 

depicted  in  Figures  12  and  13. 

The  information  in  Table  III  indicates  that  the  Jacobi  method  results  in  a 


11 


maximum  of 


-    |Vfi|l 


WZ\ 


greater  than  one;  this  value  does  not  appear  to  depend  on  grid 


position.  In  other  words,  there  is  at  least  some  component  of  the  error  that  is  magni- 
fied, instead  of  reduced,  by  this  iteration  scheme.  Moreover,  it  appears  that  the  value 
increases  with  the  number  of  intervals  on  the  grid.  Additionally,  Figures  12  and  13 
indicate  that  the  maximum  of  this  ratio  occurs  for  both  low  and  high  frequencies. 
riius.  it  appears  that  the  Jacobi  method  will  not  be  effective  in  generating  a  solution 
lu  the  point-source  problem. 


Low  frequency 

Mixed 
Frequency 

Low  frequency 

Mixed 
Frequency 

HJgh  Frequency 

Mixed 
Frequency 

Low  frequency 

Mixed 
Frequency 

1 

Low  frequency 

Figure  11.  A  two-dimensional  representation  of  frequencies  corresponding  to  a  single 
set  of  non-centered  sample  points. 
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In  order  to  compute  the  ratio     M-j    for  the  Gauss-Seidel  relaxation  scheme,  we 


rewrite  Equation  IV. 23  as 


h        n    (new)  CiKfl         (new) 


I, ml 


(new) 


(new) 


(new) 


g-m 
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'k,j 
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1.2-1 


0.8- 


0.6- 


0.4- 


J2H 


35  30 


0       0 


Figure  12.    Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  Jacobi  relaxation: 
.V  =  32,  a  =  |,  j '  =  yV  —  1.  Frequency  ranges  are  as  in  Figure  11. 


w 


here  .4  =  16J-5,  B  =  32j+5,  C  =  32j  — 11,  D  =  224j,  £  =  32.7  +  11,  F  =  32j-5,  G  = 

16j  +  5,  and  H  =  2j,  I  =  —I  +  2j,J  =  —$j,K  =  1  +  2j,  L  =  2j.  Expanding  in  a 
DFT.  equating  terms  in  the  (double)  sum  by  virtue  of  the  orthogonality  property  of 
I  lie  complex  exponential  (see  Appendix  A),  and  dividing  by  C(l,m),  we  have 
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or,  regrouping  we  have 
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1.2-1 
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0.2- 


Figure  13.    Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  Jacobi  relaxation: 
.V  =  32,  a  =  I,  j  =  N  —  1.   Frequency  ranges  are  as  in  Figure  11. 
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Thus,  the  ratio  -^  is  given  by 
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■  2ti-I 
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(IV.26) 


As  above,  we  now  seek  to  maximize 


|V/o)| 

1     i,ml 


over  the  values  I,  m  =  —  y  +  1  :  y  in  order 
to  determine  a  bound  on  how  well  this  scheme  performs.  We  again  use  Matlab  to 
compute  the  maximum  of  this  ratio  for  /,  m  =  0  :  N  for  N  =  16,32;  a  =  |,1;  and 

j  =  1,  y,  iV  —  1;  the  results  are  indicated  in  Table  IV.  Additionally,  the  matrices  of 

\v{n)\ 
values  of     '^    for  j  =  /V  —  1  are  depicted  in  Figures  14  and  15. 

The  information  in  Table  IV  indicates  that  the  Gauss-Seidel  method  results 

in  a  maximum  of     '^    less  than  one.    Thus,  all  frequency  components  of  the  error 

I     l,m\ 

are  reduced  by  this  relaxation  scheme,  which  is  what  we  seek.    However,  the  value 
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N  =  16 

N  =  32 

J  =   1 

\    \  ■ 

j=   1 

J  =  4 

J  =  /V  -  1 

I 

.8986 

.  S796  1   .8781 

.9152 

.8977 

.8970 

\  a  =   1 

.9478 

.9374    .9366 

.9567 

.9473    .9469 

H/,n)| 
fable  IV.  Maximum  Values  of     ':"1    for  /.  m  =  0  :  Ar  for  the  Gauss-Seidel  Method  oi 
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increases  with  grid  size  N,  and  also  as  the  time  stepping  is  shifted  from  Crank- 
Nicholson  to  fully  implicit.  Additionally,  while  all  of  the  values  are  less  than  one. 
as  we  move  to  grids  with  larger  N  and  toward  a  fully  implicit  time  scheme,  they 
become  dose  to  unity.  This  means  that,  although  we  may  reasonably  expect  to 
see  this  relaxation  process  converge,  it  may  be  quite  slow.  Another  point  is  that 
the  maximum  value  depends  on  grid  position;  the  shorter  the  radial  distance,  the 
larger  the  maximum,  which  apparently  is  a  reflection  of  the  radial  bias  of  the  stencils. 
Additionally,  Figures  14  and  15  indicate  that  the  maximum  of  this  ratio  occurs  only 
over  the  low  frequencies,  and  that  the  values  associated  with  error  components  over 
the  high  frequencies  are  quite  low.  This  performance  over  the  high  frequencies  will 
be  of  importance  in  Chapter  VI.  Thus,  it  appears  that  the  Gauss-Seidel  method  will 
lx1  effective  in  generating  an  iterative  solution  to  the  point-source  problem,  albeit 
potentially  slow  to  converge. 

We  now  apply  similar  techniques  to  the  radial/vertical  line  relaxation  schemes. 
The  error  equations  come  from  rewriting  Equations  IV. 20  and  IV. 21  respectively  as 
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Figure  14.   Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  Gauss-Seidel  relax- 
ation:  .V  =  32,  a  =  -,  j '  =  N  —  1.   Frequency  ranges  are  as  in  Figure  11. 
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As  above,  we  expand  these  relations  in  a  DFT,  make  use  of  the  orthogonality  property 
to  equate  terms  in  the  sums,  and  divide  by  C{l,m)  to  obtain 
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Figure  15.   Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  Gauss-Seidel  relax- 
ation: .V  =  32,  a  =  1,  j  =  N  —  1.  Frequency  ranges  are  as  in  Figure  11. 
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N=  16 

N  =  32 

i=  i 
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Table  V.  Maximum  Values  of     ' '"    for  /,m  =  0  :  A'r  for  Radial  Line  Relaxation. 
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As  above,  we  now  seek  to  maximize 


[He~n~  +  Ie     n    +J+Lc]. 


(n), 


over 


N 


the  values  /,m  =  -f  +  1  :  f.   We 


continue  to  use  Matlab  to  compute  the  maximum  of  these  ratios  for  l,m  =  0  :  N  for 
N  =  16,32;  a  =  0,  \,  1;  and  j  =  1,  y,  N  —  I;  the  results  are  indicated  in  Tables  V 


M 


Hi 


and  VI.   Additionally,  the  matrices  of  values  of     'ffi    for  j  =  N  —  I  are  depicted  in 
Figures  16,  17,  18,  and  19. 

The  information  in  Tables  V  and  VI  indicates  that  both  line  relaxation  meth- 

\v{n)\ 
ods  result  in  a  maximum  of     '{^    less  than  one;  radial  line  relaxation  seems  to  promise 
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for  /,  m  =  0  :  Ar  for  Vertical  Line  Relaxation. 


better  performance.  Many  of  the  results  for  the  line  relaxation  schemes  parallel  those 
noted  for  the  Gauss-Seidel  method.  All  frequency  components  of  the  error  are  reduced 
by  the  line  relaxation  schemes;  their  performance  would  appear  to  be  slightly  better 
than  i  hat  for  the  Gauss-Seidel  method.  Additionally,  the  maximum  values  increase 
with  both  the  grid  size  N.  and  as  the  time  stepping  is  shifted  from  Crank-Nicholson 
to  fully  implicit.  As  with  the  Gauss-Seidel  method,  while  all  of  the  values  are  less 
than  one,  as  we  move  to  grids  with  larger  N  and  toward  a  fully  implicit  time  scheme, 
they  become  close  to  unity.  This  means  that,  although  we  may  reasonably  expect 
to  see  this  relaxation  process  converge,  it  may  be  quite  slow,  and  that  the  increased 
amount  of  computational  work  as  compared  with  work  for  the  Gauss-Seidel  method 
may  not  be  worth  the  payoff.  The  maximum  values  for  these  schemes  also  depend  on 
grid  position.  For  the  horizontal  line  relaxation,  the  shorter  the  radial  distance,  the 
smaller  the  maximum;  for  the  vertical  line  relaxation,  the  shorter  the  radial  distance, 
the  larger  the  maximum.  Once  again,  apparently,  this  is  a  reflection  of  the  radial  bias 
of  the  stencils.  Additionally,  Figures  16,  17,  18,  and  19  indicate  that  the  maximum 
of  these  ratios  occurs  only  over  low  frequencies,  and  that  the  values  associated  with 
error  components  over  the  high  frequencies  are  quite  low.  In  fact,  the  values  for  the 
line  relaxation  schemes  over  the  high  frequencies  appear  to  be  about  one-half  the  cor- 
responding values  for  the  Gauss-Seidel  method.  As  noted  earlier,  this  performance 
over  the  high  frequencies  will  be  of  importance  in  Chapter  VI. 

Some  general  trends  are  evident  from  the  results  of  the  numerical  experiments. 
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Figure  16.  Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  radial  line  relaxation: 
N  =  -Y2,  a  =  ^,  j  =  N  —  1.  Frequency  ranges  are  as  in  Figure  11. 

Since  the  maxima  should  all  be  less  than  one  to  guarantee  convergence,  it  appears  that 
the  Jacobi  method  will  not  be  useful  in  iteratively  solving  the  point-source  problem 
(Equation  II. 2).  Of  the  schemes  that  have  maxima  less  than  one,  those  that  have  the 
smallest  maxima  are  those  that  reduce  error  components  most  quickly  and,  therefore, 
arc  the  schemes  that  will  converge  most  quickly.  The  line  relaxation  schemes  have  the 
smallest  maxima;  the  maximum  values  over  the  high  frequencies  are  about  one-half 
the  corresponding  values  for  the  Gauss-Seidel  method,  but  there  is  little  increase  in 
performance  over  the  low  frequencies.  To  determine  whether  or  the  extra  work  in 
using  line  relaxation  is  worth  the  effort,  a  cost  comparison  of  getting  to  the  same 
level  of  error  as  other  relaxation  schemes  will  need  to  be  done.  Additionally,  the 
maximum  values  vary  depending  on  the  time  stepping  scheme  and  the  grid  position; 
they  are  higher  for  implicit  time-stepping  than  for  Crank-Nicholson,  indicating  that 
convergence  will  take  longer  for  the  former  as  compared  to  the  latter,  and  the  maxi- 
mum values  increase  with  N.  Hence,  as  we  move  to  a  grid  with  more  nodes,  we  not 
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Figure  17.   Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  radial  line  relaxation: 
V  =  32,  a  =  1,  j  =  N  —  1.  Frequency  ranges  are  as  in  Figure  11. 

only  have  more  work  to  do  by  virtue  of  the  greater  number  of  equations  to  solve,  we 
also  have  to  work  harder  to  reduce  each  component  of  the  error.  That  is,  the  effect 
of  moving  to  a  finer  grid  is  to  more  than  quadruple  the  work  required  to  solve  the 
problem. 

The  Gauss-Seidel  and  line  relaxation  schemes  appear  to  be  the  most  efficient 
in  terms  of  reducing  error  components.  Additionally,  in  Chapter  VI,  the  performance 
of  these  schemes  over  the  high  frequencies  will  be  of  interest;  we  will  want  to  use  those 
schemes  that  have  the  smallest  maximum  value  over  the  high  frequencies.  In  particu- 
lar, for  the  Gauss-Seidel  and  line  relaxation  schemes,  the  high  frequency  components 
of  the  error  are  eliminated  much  more  quickly  than  the  low  frequency  components. 
Considering  the  amount  of  work  to  be  done  and  the  performance  of  the  scheme,  we 
implement  the  Gauss-Seidel  method  in  the  solution  process. 
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Figure  IS.   Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  vertical  line  relax- 
ation: i'V  =  32,  a  =  i  j '  =  N  —  1.  Frequency  ranges  are  as  in  Figure  11. 


o    o 


Figure  19.   Ratio  of  amplitudes  of  new  to  old  Fourier  coefficients,  vertical  line  relax- 
ation: /V  =  32,  a  =  1,  j :  =  N  —  1.  Frequency  ranges  are  as  in  Figure  11. 
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V.  ITERATIVE  SOLUTION 

The  results  of  Chapter  IV  indicate  that  the  Gauss-Seidel  method  is  the  sim- 
plest of  the  relaxation  schemes  that  appear  to  result  in  a  converging  sequence  of 
approximate  solutions,  and  therefore  Gauss-Seidel  is  our  method  of  choice.  Having 
chosen  a  relaxation  scheme,  we  must  now,  as  part  of  the  iterative  solution  process, 
make  suitable  choices  for  the  type  of  time-stepping,  a,  and  the  time  step  size,  g\  the 
goal  is  to  maximize  accuracy  and  minimize  computational  work.  We  now  conduct 
experiments  to  solve  Equation  IV. 3  by  iterating  with  a  solver  which  is  based  on  the 
iteration  matrix.  Pry  —  —  (D  -f  L)-1U,  and  use  this  process  to  evaluate  some  of  the 
possible  choices  for  a  and  g.  The  technique  is  to  solve  first  at  a  single  time  step  by 
iterating  until  the  difference  between  successive  approximations,  ||77J+l»  —  Tl\\\  or, 
perhaps  better,  ||MT^+1)  —  /(T/M||,  is  negligible.  This  approximate  solution.  Tn,  is 
used  as  input  for  the  time-stepping  scheme  so  that 

{N+V2 


f(Tn)  =   £  [-  /  WdV  +  (!-«)«/  v<#1  •  hdS]T; 

,.      a  Jv  Js 


1 
i=i     9 

(as  in  Chapter  III)  becomes  the  new  right  hand  side  in  Equation  IV. 3.  The  iteration 
process  is  repeated  to  obtain  the  approximation  at  the  new  time  step.  The  solution 
is  stepped  in  time  in  this  fashion  until  the  difference  between  solutions  at  successive 
time  steps  is  negligible,  representing  a  steady  solution. 

An  algorithm  for  this  process  might  look  something  like  the  following: 
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Given  an  initial  guess.  T°.  and  tolerances  Si,  S2, 

Set  n  =  0 

While  ||fn+1  -  fn\\  >  S2 

1)  Compute  f{fn) 

2)  Set  .s  =  1,  fjj1  =  fn 

3)  While  Hf^1-- 77^1,11  ><Ji 

77;+;,  f-  Pcf,^1  4-  (D  +  L)-1/!^) 

.s  =  5  +  l 
End  while 

4)  fn+i   <_  fn+i 

Return  fn+1  as  the  solution  at  t  =  (n  +  I) At 
n  =  n  +  1 
End  while 

^steady    , "T>n 

Stop 


Starting  with  an  initial  temperature  distribution  of  T  —  0  everywhere,  bound- 
ary conditions  specified  by  Equations  11.18,  11.19,  and  11.20,  iV  =  16,  and  a  =  ^, 
several  values  of  g  are  used  in  an  attempt  to  achieve  a  reasonable  solution.  The  re- 
sults of  these  experiments  indicate  that  for  larger  values  of  g,  say  g  =  h,  the  solution 
exhibits  a  non-physical  oscillation  in  time  (see  Figure  20)  at  the  origin,  where  the 
heat  source  is  located.  Instead,  the  solution  should  monotonically  increase  until  it 
levels  olf  at  steady  state. 

The  oscillation  can  be  removed  by  making  the  temporal  step  size  smaller,  say 
g  =  /?'3  (see  Figure  21).  This  step  size,  however,  results  in  negative  temperatures 
along  the  z  =  0  axis  near  the  origin.  Negative  temperatures,  like  the  oscillations,  are 
physically  meaningless.  A  valid  solution  process  must  eliminate  both  of  these  non- 
physical  characteristics  from  the  solution.  However,  it  turns  out  that  for  a  =  t,  there 
is  no  value  of  g  that  produces  a  physically  meaningful  solution.  In  particular,  for 
g  =  /)2,  there  are  both  oscillations  and  negative  values  in  the  approximate  solution 
(see  Figure  22  for  a  comparison  of  several  values  of  g).  Thus,  in  order  to  compute  a 
realistic  solution,  it  is  necessary  to  use  values  of  a  >  |. 
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Maximum  Solution  Values 


Figure  20.  The  maximum  values  of  the  solution  for  each  time  step,  where  a  =  -  and 
g  =  h. 

In  an  attempt  to  produce  a  solution  that  is  physically  meaningful,  numerical 
experiments  are  conducted  using  values  of  ^  <  a  <  1,  and  values  of  g  ranging  from 
/r  to  h,  for  each  value  of  a.  Solutions  with  no  oscillation  and  no  negative  values  are 
possible  for,  say  a  =  |,  but  this  requires  that  the  temporal  step  size  be  as  small  as 
g  =  h~sTT .  The  size  of  the  temporal  step  indicates  how  much  work  will  have  to  be 
done  to  reach  a  steady  solution;  the  smaller  the  step,  the  longer  to  reach  steady  state. 
In  order  to  minimize  the  work,  we  must  maximize  the  temporal  step  size.  Further 
experimentation  indicates  that  for  a  =  1  (fully  implicit  time  stepping)  and  g  =  h, 
the  result  is  not  only  a  physically  meaningful  solution  (i.e.,  no  oscillations  and  no 
negative  values),  but  also  a  reasonable  temporal  step  size.  Therefore,  the  iterative 
solution  is  computed  using  these  values. 

Iterative  solutions  for  N  =  8, 16,  and  32,  are  computed  using  the  same  initial 
temperature  distribution  and  boundary  conditions  as  above.  One  measure  of  accu- 
racy of  these  solutions  is  determined  by  comparing  them  with  an  analytic  solution 
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Maximum  Solution  Values 


Figure  21.  The  maximum  values  of  the  solution  for  each  time  step,  where  a  =  -  and 
g  =  h\ 

(Equations  11.15  and  11.16).  However,  the  problem  that  gives  rise  to  the  analytic  so- 
lution presented  in  Chapter  II  is  somewhat  different  from  the  linear  algebra  problem 
that  we  are  solving.  Therefore,  as  in  [Ref.  11]  we  approximate  the  exact  solution  of 
Equation  II. 2  by  solving  Equation  IV. 3  for  N  =  64  (see  Figure  23).  We  then  compare 
this  solution  with  solutions  of  Equation  IV. 3  on  coarser  grids  (N  =  32  and  N  =  16) 
to  get  approximations  of  the  discretization  errors,  which  may  be  used  to  give  an  in- 
dication of  the  order  of  the  accuracy  of  the  solution.  The  measure  that  we  use  is  the 
discrete  energy  norm,  defined  by 


|Zr|||  =  {LhD\Dh)*, 


(V.l) 


where  (•,•}  denotes  the  Euclidean  inner  product,  Lh  is  the  operator  defined  in  Equa- 
tion IV. 2,  and  Dh  =  Th  —  T*  approximates  the  discretization  error,  where  Th  is  the 
solution  to  Equation  IV. 2  on  grid  k,  and  T*  is  the  "exact"  solution  (i.e.,  solution  of 
Equation  IV. 2  with  N  =  64  sampled  on  grid  h).  (see  [Ref.  11]) 
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Maximum  Solution  Values 


8  10  12 

Time  steps 

Figure  22.  The  maximum  values  of  the  solution  for  each  time  step,  where  a  —  |;  uo" 
indicates  g  =  /i,  "+"  indicates  g  —  h2,  "*"  indicates  g  =  /12,  "-"  indicates  g  =  /i3. 
Care  must  he  taken  in  analyzing  the  figure  since  the  apparent  equivalence  of  time 
steps  can  be  misleading,  e.g.,  for  g  =  h2  and  h  =  j?,  it  takes  16  time  steps  to  ''equal" 
one  time  step  for  g  =  h. 


The  solutions  for  an  inital  time  step1,  for  a  time  of  one  second,  and  for  steady 
state  have  been  compared  using  the  discrete  energy  norm,  with  the  resulting  dis- 
cretization errors  indicated  in  Table  VII.  We  might  hope  to  see  a  common  factor  of 
decrease  as  we  move  to  finer  grids.  That  is,  if  the  discretization  error  on  grid  h  were 
of  the  order  e  =  chp,  for  some  constant  c,  then  the  ratio  of  errors  for  grids  h  and  ^ 


1  With  y  —  h,  the  size  of  a  time  step  differs  with  grid  size.  We  have  adjusted  for  this  difference  by 
requiring  that  the  same  "amount"  of  time  elapse  for  solutions  used  in  the  comparison,  e.g.,  for  the 
initial  time  step,  the  "exact"  solution  on  N  =  64  after  eight  time  steps  is  compared  to  the  solution 
on  N  =  32  after  four  time  steps,  to  the  solution  on  N  =  16  after  two  time  steps,  and  to  the  solution 
on  N  =  8  after  one  time  step. 
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Steady  Solution  for  N=64 


0     0 


Figure  23.  The  steady  state  "exact"  solution  for  N  =  64. 
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Thus,  a  decrease  of  a  factor  of  four  would  indicate  that  the  process  achieves  0(h2) 
discretization  error  (see  [Ref.  11]).  Table  VII  indicates  that  the  discretization  error 
decreases  by  a  factor  of  about  two  when  moving  from  grid  size  N  =  8  to  N  =  16, 
and  by  a  factor  of  about  four  when  moving  from  grid  N  =  16  to  N  =  32.  A  grid 
size  N  —  8  may  be  too  small  to  adequately  reflect  the  accuracy  of  the  solution,  which 
would  leave  the  factor  of  four,  possibly  suggesting  that  the  discretization  error  is 
0(li2).  Even  so,  a  specific  estimate  of  the  discretization  error  is  best  deferred  until 
more  information  can  be  obtained,  i.e.,  solutions  on  finer  grids. 
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N  =  8 

.'V  =  16 

N  =  32 

Initial  time  step 

43.07 

24.02 

6.85 

Time  =  "1  second" 

42.99 

23.98 

6.87 

Steady  state 

42.99 

23.98 

6.87 

Table  VII.  Discretization  Errors  for  the  Iterative  Solution. 
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VI.         MULTILEVEL  SOLUTION 

Now  that  we  can  solve  a  discrete  representation  of  the  heat  equation  (Equa- 
tion IV. 2)  iteratively,  we  want  to  apply  multigrid  techniques  in  order  to  attempt  to 
accelerate  the  solution  process.  Multigrid  is  a  method  to  improve  on  solution  by 
relaxation  by  making  use  of  the  advantages  of  working  on  successively  coarser  grids 
(for  a  more  complete  treatment,  see  [Ref.  3]).  It  has  been  used  with  marked  suc- 
cess in  solving  a  variety  of  problems,  specifically  elliptic  partial  differential  equations. 
Although  the  problem  we  are  solving  is  parabolic,  in  the  time-stepping  regime  we 
must  solve  an  elliptic  problem  at  each  time  step.  Thus,  multigrid  may  prove  useful 
in  streamlining  our  solution  process.  We  begin  by  introducing  the  multigrid  method 
and  some  of  the  exigencies  of  the  use  of  multiple  grids.  We  then  analyze  some  of  the 
characteristics  of  the  method  using  local  mode  analysis,  and  present  numerical  re- 
sults from  the  solution  process.  The  amount  of  work  to  achieve  the  iterative  solution 
serves  as  a  baseline  against  which  the  computational  work  for  the  multilevel  solution 
is  compared. 

As  indicated  in  the  analysis  of  the  Gauss-Seidel  and  line  relaxation  schemes  in 
Chapter  IV,  the  high  frequency  (oscillatory)  components  of  the  error  are  extirpated 
much  more  efficiently  by  these  schemes  than  are  the  low  frequency  (smooth)  compo- 
nents. The  result  of  applying  a  relaxation  scheme  to  generate  an  approximate  solution 
on  a  specific  grid  size  (call  it  a  fine  grid,  size  N),  is  that  the  oscillatory  components 
of  the  error  are  smoothed.  After  sufficient  work  has  been  done  on  the  fine  grid  to 
smooth  the  error  (in  effect,  when  the  relaxation  process  stalls),  the  problem  may  be 
shifted  to  a  coarser  grid  (grid  size  y)  where  the  smooth  components  of  the  error 
become  oscillatory  (see  figure  24,  taken  from  [Ref.  ',]]).  The  relaxation  scheme  is  then 
applied  again  to  smooth  the  oscillatory  components  of  the  error.  The  information 
gained  from  smoothing  on  the  coarse  grid  is  transferred  back  to  the  fine  grid,  where  it 
becomes  a  correction  to  the  original  approximation.  That  is,  the  result  of  smoothing 
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mi  the  fine  grid,  transferring  to  the  coarse  grid,  smoothing  on  the  coarse  grid,  and 
transferring  back  to  the  fine  grid  is  to  produce  a  coarse  grid  correction.  This  pro- 
cess of  using  multiple  grids  to  obtain  an  approximate  solution  has  the  advantages  of 
a)  more  effectively  targeting  the  error  components,  and  b)  requiring  less  work,  since 
the  coarse  grid  has  only  2  the  number  of  points  on  the  fine  grid,  where  d  is  the 
dimension  of  the  domain.  Multigrid  also  allows  for  ways  to  improve  upon  the  initial 
guess  that  is  used  in  the  smoothing  (see  [Ref.  '■)]).  However,  since  we  are  solving  a 
time-stepping  problem,  where  the  current  solution  becomes  the  initial  guess  for  the 
next  time  step,  this  will  not  be  a  concern  for  us.  Moreover,  while  there  is  a  good  deal 
more  to  multigrid  than  coarse  grid  correction,  this  concept  forms  the  foundation  of 
our  solution  process. 
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k  =  4  wave  on  N  =  16 


a 


k  =  4  wave  on  N  =  8 


Figure  24.  A  wave  with  wave  number  k  =  4  on  Qh  (N  =  16)  is  projected  onto  Q,H 
(N  =  8).  For  N  =  16,  k  —  4  implies  that  the  wave  is  -^  =  ^  the  way  up  the  spectrum, 
while,  for  N  =  8,  k  =  4  is  the  wave  that  is  |  =  |  way  up  the  spectrum.  Thus,  the 
coarse  grid  "sees"  a  wave  which  is  more  oscillatory. 

In  order  to  make  effective  use  of  this  technique,  we  must  specify  what  informa- 
tion to  transfer,  and  how  to  transfer  it;  first  we  consider  what  information  to  transfer. 
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R.erall  that  we  are  attempting  to  solve  Equation  IV. 3. 

Ml  =  /, 

.mil  suppose  that  we  have  an  approximate  solution  v\  the  error  is  given  bye=  T*  -v, 
where  /'*  is  the  exact  solution.  Thus,  the  error  satisfies 

Mc    =    / -  Mv 

=    f,  (VI.  1) 

where  r  is  the  residual.  The  smooth  components  of  the  error  are  the  troublesome 
ones,  since  the  relaxation  schemes  "kill"  the  oscillatory  modes.  Since  smooth  modes 
on  <i  line  grid  become  oscillatory  when  projected  onto  a  coarse  grid,  it  is  natural  to 
consider  transferring  a  representation  of  the  error,  i.e.,  the  residual,  from  one  grid  to 
the  next1.  In  this  way,  the  relaxation  schemes  used  on  the  coarse  grid  will  reduce 
components  of  the  error  that  could  not  be  reduced  on  the  fine  grid.  Additionally,  we 
know  that  relaxation  smooths  the  error,  and  therefore  we  can  accurately  represent 
the  error  on  the  coarse  grid.  After  transferring  the  residual  to  the  coarse  grid,  we 
can  relax  on  the  coarse  grid  version  of  the  residual  equation  (Equation  VI. 1),  solving 
for  the  error.  Thus,  when  we  have  determined  an  approximation  to  the  error  on  the 
coarse  grid,  we  can  transfer  it  back  to  the  fine  grid  as  a  correction  to  the  current 
approximation.  That  is,  since  T*  =  v  -f  e,  if  we  know  v  we  can  correct  it  by  adding 
an  approximation  of  e.  This  process  can  be  outlined  as  in  the  following  steps,  where 
h  represents  the  grid  spacing  on  the  fine  grid,  H  the  spacing  on  the  coarse  grid,  and 
we  assume  2h  =  H . 

Coarse  grid  correction  (two-grid  scheme) 

•   Relax  on  M/lT/l  =  fh  on  iVl  to  obtain  an  approximation  vh. 


'There  are  a  number  of  other  good  reasons  to  transfer  the  residual,  as  outlined  in  [Ref.  3]. 
However,  there  are  other  multigrid  schemes  that  do  not  transfer  the  residual,  such  as  the  Full 
Approximation  Scheme  (FAS)  (see  [Ref.  2]  and  [Ref.  1]).  FAS  is  useful  in  dealing  with  nonlinear 
problems  but,  since  our  problem  is  linear,  we  will  not  employ  it. 


(i.-i 


Compute  the  residual  r*  =  fh  -  M 


htfh 


•  Solve  the  residual  equation  Mr      —  r      on  fr     to  obtain  an  approximation 
to  the  error  e    . 

•  Correct  the  approximation  obtained  on  Q,h  with  the  error  estimate  obtained 
on  ilH:  vh  f-  u^  +  e^Ref.  3] 

The  task  of  solving  on  Q,h  can,  thus,  be  exchanged  for  the  task  of  relaxing  on 
Hk.  and  then  solving  on  Q,H .  The  requirement  to  solve  on  QH  can  be  treated  in  like 
fashion;  the  task  of  solving  can  be  exchanged  for  the  task  of  relaxing,  and  then  solving 
on  the  next  coarser  grid.  This  procdeure  of  transferring  to  successively  coarser  grids 
can  be  continued  until  a  grid  is  reached  on  which  an  "exact"  solution  is  possible. 
That  is,  the  solution  can  be  generated  by  recursive  use  of  the  coarse  grid  correction, 
which  can  be  represented  as  follows: 

Solve  on  Q,'1  by  relaxing  on  tth  and  solving  on  02/l. 

Solve  on  Vt2k  by  relaxing  on  il2k  and  solving  on  f24/l. 

Solve  on  Q,4h  by  relaxing  on  £74/l  and  solving  on  Q,8  . 

"Exact"  solve  on  coarsest  grid. 

Correct  on  ft4'1. 
Correct  on  H'2/l. 
Correct  on  $Vl. 

This  process  of  starting  on  a  fine  grid,  moving  to  the  coarsest  grid,  and  then  going 
back  to  the  fine  grid  is  known  as  a  V-cycle  (see  Figure  25).  Each  time  the  transfer 
to  a  coarser  grid  is  made,  the  relaxation  process  there  smooths  another  portion  of 
the  error  component  spectrum  (see  Figure  26).  Thus,  by  the  time  the  coarsest  grid 
is  reached  and  the  problem  has  been  solved,  all  of  the  error  components  have  been 
smoothed.  When  the  problem  is  solved  on  the  coarsest  grid,  the  approximate  error 
is  transferred  to  the  next  finer  grid  to  become  a  correction.  The  corrected  error  is 
then  passed  to  the  next  finer  grid,  and  so  on  until  we  are  back  on  the  finest  grid, 
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where  the  original  approximation  is  updated  with  the  composite  correction  from  the 
i  ii.ii sci  grids.  While  there  are  other  multigrid  methods  which  may  prove  more  useful, 
depending  on  the  nature  oi  the  problem  ( see  [Ref.  3]),  we  use  the  V-cycle. 

h 

2h 
4h 
8h 
16h 

Figure  25.  Schedule  of  grids  for  a  V-cycle. 

HhH 1 1 1 1 

h      mmm—^m^——^—^^^^—^^^^— 


ih 


4h 


Xh 


HM 1 1 1 1 

Figure  26.  Error  component  spectrum  for  various  grid  sizes.  On  the  Hue  grid,  //. 
relaxation  smooths  the  high  frequency  components  of  the  error  (heavy  black  line). 
When  the  problem  is  transferred  to  the  next  coarser  grid,  2/i,  the  high  frequencies 
of  the  remaining  (unsmoothed)  components  (heavy  black  line)  are  smoothed.  This 
process  continues  all  the  way  down  to  the  coarsest  grid,  where  the  problem  is  solved, 
eliminating  the  remaining  components  of  the  error.  By  continuing  to  transfer  to 
successively  coarser  grids,  all  frequency  components  of  the  error  are  smoothed. 

So  far,  we  have  not  indicated  how  information  will  transferred,  nor  how  an 
error  estimate  on  the  coarse  grid  (solving  MH cH  =  fH)  will  be  computed.    We  will 
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deal  with  the  former  question  in  the  next  section;  the  latter  question  revolves  around 
deciding  what  to  use  as  the  coarse  grid  operator,  M  .  One  approach  is  to  simply 
discretize  the  problem  on  the  coarse  grid  in  precisely  the  same  fashion  that  it  is 
discretized  on  the  fine  grid.  This  has  the  advantages  that  the  work  has  already  been 
dune,  and  that  it  is  easy  to  implement  computationally.  While  this  method  generally 
works  (on  model  problems),  it  has  the  disadvantage  that  it  is  often  not  true  to  the 
physics,  e.g.,  conservation  may  no  longer  be  enforced.  There  are  other  methods  for 
determining  the  coarse  grid  operator  which  are  true  to  the  physics  and  allow  for  strong 
theoretical  treatment  of  convergence  and  other  properties  (see  [Ref.  2]  and  [Ref.  3]). 
However,  because  discretization  on  the  coarse  grid  generally  works,  and  because  it 
simplifies  the  computations,  forming  Mw  in  this  manner  is  the  method  that  we  use. 
We  conclude  this  section  with  a  discussion  of  what  happens  to  boundary  con- 
ditions when  we  transfer  the  residual.  The  boundary  conditions  for  our  problem  are 
a  mixture  of  Dirichlet  and  Neumann  conditions  (Equations  11.18,  11.19,  and  11.20). 
By  relaxing  on  Equation  IV. 2, 

L\T]  =  /, 
(where  we  have  dropped  the  superscript  h)  we  obtain  an  approximation  u,  which  gives 
rise  to  the  residual  equation,  r  =  f  —  L[v].  Suppose  that  T*  is  the  exact  solution,  so 
that  /  =  L[T*]\  we  require  that  the  approximation  satisfy  the  same  conditions  as  the 
exact  solution  on  the  boundary,  dtt.  Thus,  the  residual  equation  becomes 

r    =     L[T*)  -  L[v] 

=     [  f  (-  -  a/cV2)rW]  -[/(--  aKV2)vdV] 
Jv  g  Jv  g 

=     -  [  (T*  -  v)dV  -ok  I  (V2T*  -  V2v)dV 
g  Jv  Jv 

=    -  [  (T*  -  v)dV  -an  f  (VT*  -  Vi;)  •  MS.  (VI.2) 

g  Jv  Js 

Therefore,  since  the  boundary  conditions  satisfy  (T*  =  v)\dn  and  (^-  =  g^)|ao2,  all 


2 -$-  —  V  •  h  denotes  the  outward  normal  on  dQ. 

on 
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boundary  conditions  for  7'  become  homogeneous;  the  first  term  on  the  right  hand  side 
ol   Equation  VI. 2  indicates  that  this  is  true  tor  the  Dirichlet  conditions,  the  second 

term  indicates  that  it  is  true  tor  the  Neumann  conditions.  Thus,  since  we  are  solving 
the  residual  equation  on  all  the  successively  coarser  grids,  we  impose  homogeneous 
boundary  conditions  on  all  but  the  finest  grid. 

A.      INTERGRID  TRANSFERS 

In  order  to  transfer  information  between  grids,  we  develop  operators  that 
restrict  data  from  a  tine  grid  to  a,  coarse1  grid,  and  interpolate  data  from  a  coarse 
grid  to  a  fine  grid  (see  [Ref.  '.]]).  These  operators  are  designated  by  ^"bsci-ipt1"-  wnere 
the  subscript  indicates  the  grid  from  which  the  information  is  being  transferred,  and 
the  superscript  indicates  the  grid  to  which  the  information  is  being  transferred.  Thus, 
i  he  restriction  operator  is  indicated  by  I//,  since  the  information  goes  from  the  line 
grid  to  the  coarse  grid,  and  the  interpolation  operator  is  represented  by  1^,  since 
the  information  goes  from  the  coarse  grid  to  the  fine  grid.  In  other  words,  restriction 
is  a  process  of  determining  what  values  to  assign  to  coarse  grid  points,  based  on  the 
values  at  fine  grid  points;  interpolation  is  a  process  of  determining  what  values  to 
■  issign  to  fine  grid  points,  based  on  the  values  at  coarse  grid  points.  In  [Ref.  1],  such 
operators  are  developed  which  take  advantage  of  FVE  characteristics.  We  consider 
two  types  of  restriction  operators,  one  which  follows  from  the  physics  and  one  which 
is  very  simple,  and  one  type  of  interpolation  operator. 

One  of  the  advantages  of  the  FVE  method  is  its  fidelity  to  conservation  consid- 
erations. The  restriction  technique  that  we  describe  first  follows  from  the  conservation 
notion.  In  order  to  determine  what  value  to  assign  to  a  coarse  grid  point,  the  fine  grid 
is  laid  over  the  coarse  grid,  where  the  coarse  grid  control  volumes  align  with  fine  grid 
lines  (see  Figure  27).  The  control  volume  for  a  given  coarse  grid  point  includes  all  or 
part  of  the  control  volumes  of  nine  fine  grid  points.  Hence,  the  value  of  a  quantity  on 
the  coarse  grid  includes  contributions  from  the  value  of  the  quantity  on  each  of  the 
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nine  fine  grid  points.  For  each  of  these  fine  grid  points,  the  contribution  is  determined 
by  computing  the  fraction  of  the  associated  fine  grid  volume  that  is  contained  in  the 
coarse  grid  volume.  The  value  of  the  quantity  at  the  fine  grid  point  is  multiplied  by 
this  fraction,  and  the  result  is  the  contribution  to  the  coarse-grid  value.  The  value 
of  the  quantity  at  the  coarse  grid  point  is  thus  the  sum  of  the  nine  contributions 
determined  in  this  manner. 

Fine  Grid 


grid  lines 


control  volumes 


Coarse  Grid 


Figure  27.  The  "conservation  restriction"  operator. 

In  the  case  of  cylindrical  coordinates,  care  must  be  taken  when  deciding  what 
percentage  of  the  fine  grid  volumes  fall  within  each  coarse  grid  volume.  Since  the 
volumes  change  with  the  radius,  the  restriction  operator  will  be  a  function  of  grid 
location.  The  stencil  for  the  restriction  operator  gives  the  percentage  of  the  fine  grid 
value  that  is  assigned  to  the  coarse  grid  value  based  on  the  percentage  of  the  fine 
grid  volume  that  falls  within  the  coarse  grid  volume  (see  Figure  28).  For  example, 
the  amount  of  the  value  of  the  point  at  (k,  j  +  1)  that  is  used  in  the  sum  is  found  by 


68 


: 

lu     i           cu          t     ru 

Ic     |           cc           '     re 

Id     !           cd          !     rd 

grid  lines 


control  volumes 


(j- nil 


jh 


(j+l)h 


Figure  28.  The  tine  grid  region,  comprising  contributions  from  the  nine  tine  grid 
control  volumes,  that  corresponds  to  a  coarse  grid  control  volume.  The  central  fine 
grid   point   is  labeled  (k,j),   (k,j  even)  for  which  the  radial  distance  r  =  jh,  and 

i  ui  responds  to  coarse  grid  point  (-,  |). 


first  determining  how  much  of  its  fine  grid  control  volume  (right-center,  or  re)  falls 
within  the  region  that  corresponds  to  the  coarse  grid  control  volume.  This  volume 
is  then  divided  by  the  total  control  volume  for  the  point  (k,j  +  1 )  to  determine  the 
volume  percentage.  So,  if  the  fine  grid  location  is  at  the  radial  distance  r  =  jh.  then 
we  have  for  the  volume  re  (a  toroidal  prism) 

1 


7rh[(j+l)2h2-(j  +  -)2h2)}, 


or 


which  is  divided  by  the  fine  grid  control  volume  for  the  point  at  (A;,  j  +  1), 


or 


7Th3(2j  +  2), 


69 


giving 


i  +  I 


2(j  +  l) 

Similarly,  the  ratio  for  the  left-center  (lc)  portion  of  the  stencil  is 

2(i-i)' 

and  the  remainder  of  the  stencil  is  given  by, 


1 1  j- 


i      i/i+7 


4  V  j  — 1  >       2        4  V j  +  1 


1^-1 
-1 

3 


!(J—±)      1      ifiJJ 


l/J  +  t 
1 
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;L_i  ) 


i     i  /j+ 


4V j-\  '       2        4Vj+l 
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This  is  a  stencil  for  interior  points;   similar  calculations  determine  the  stencils  for 
boundary  points. 
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(j-3)h 


(j-2)h 


(j-l)h 

fine  grid  lines 

fine  grid  control  volumes 


Jh 


(j+l)h 


_  _         coarse  grid  control  volumes 

0  coarse  grid  points 

Figure  29.  Neighboring  coarse  grid  points. 

Since  this  operator  is  both  somewhat  out  of  the  ordinary  and  dependent  on  grid 
location,  we  consider  it  worthwhile  to  try  to  get  a  sense  of  how  this  restriction  operator 
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iicai s  the  quantities  being  restricted.  For  this  purpose,  we  consider  a  comparison 
between  the  values  assigned  to  neighboring  coarse  grid  points,  as  in  Figure  29,  when 
the  restriction  operator  acts  on  a  grid  with  constant  value  Tk,j  =  L.  The  value  for 
coarse  grid  point  number  1  is  given  by 

J  ~  1  J  +  1 


or 


2  +  2( * )     =    2  +  2(  — -^-^-r) 

i 

=     4  + 


2(i  +  l)(j-l 
he  corresponding  value  for  coarse  grid  point  number  2  is 

£ J-  +2  4-  - i 

J  -  3  j  -  1 


or 


2  +  2(-± J       4       =    2  +  2(± J  4 ) 

Hi-3)(i-i)  +  l  (i  -  3)(i  -  i) 


2(j-3)(j-l) 
Since  (j  +  1)  >  (j  —  3)  for  2  <  j  <  N  ~  2,  then  -4^-  <  -j-^,  which  implies 


1  1 

< 


2(i  +  l)(i-l)    "  2(i  -  3)(j  -  1) 

and  therefore 

4  4-  777-. 7T7-. 77   <  4  + 


2(j  + l)(j-l)  2(i  -  3)0'  -  1) 

Thus,  if  the  restriction  operator  acts  on  a  constant  vector,  the  values  assigned  to  the 
coarse  grid  points  decrease  as  the  radial  distance  increases. 

Does  this  make  sense  in  terms  of  the  physics  of  the  problem?  Should  a  constant 
function  remain  a  constant  after  restriction,  or  not?  In  general,  one  might  expect 
that  restriction  would  transfer  a  constant  to  a  constant.  The  conservaton  restriction 
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operator  acts  on  the  basis  of  volume  percentages  and,  thus,  considers  the  values 
assigned  to  each  control  volume  in  terms  of  the  actual  volume.  In  some  sense,  this 
restriction  operator  converts  the  values  assigned  to  a  control  volume  into  a  "density", 
and  then  transfers  this  amount.  In  the  case  of  cylindrical  coordinates,  assigning  a 
constant  value  to  each  control  volume  means  that  the  ratio  of  the  value  to  its  volume 
decreases  with  radial  distance.  Thus,  the  result  of  the  above  computation  is  no 
surprise.  Moreover,  it  seems  to  make  sense  physically  to  transfer  "densities"  in  terms 
of  enforcing  conservation.  However,  in  order  to  regain  the  fine-grid  interpretation  of 
the  information  after  it  has  been  transferred  may  require  more  work.  In  particular, 
some  compensation  may  need  to  be  made  relative  to  the  coarse  grid  volumes  (i.e., 
convert  the  "density7  back  to  the  original  type  of  information).  The  question  of  how 
or  whether  to  make  this  compensation  will  not  be  addressed  here.  Additionally,  it  is 
certainly  possible  that  the  information  to  be  transferred  does  not  make  sense  in  the 
context  of  densities.  In  that  case,  physical  intuition  may  need  to  be  suspended  while 
the  efficacy  of  the  solution  process  is  determined. 

The  other  restriction  operator  we  introduce  is  the  injection  operator.  Injection 
is  accomplished  by  simply  assigning  values  to  the  coarse  grid  points  directly  from  the 
corresponding  fine  grid  points,  vjfj  =  v2\2j  (see  [Ref-  3]).  The  injection  operator 
has  the  advantage  of  being  a  linear  operator,  whereas  the  conservation  restriction 
operator  may  not  be  linear.  Note  that,  when  the  conservation  restriction  operator  is 
applied  to  a  grid  with  constant  value  Tk,j  =  1,  the  result  is  T2k,2j  ~  4.  If  the  injection 
restriction  operator  is  applied  to  the  same  grid,  T<ik,2j  —  I-  Thus,  care  must  be  taken 
when  comparing  these  two  restriction  operators. 

In  order  to  transfer  from  the  coarse  grid  to  the  fine  grid,  we  must  choose  an 
interpolation  operator.  Since  we  are  assuming  that  the  basis  functions  for  the  FVE 
method  are  linear,  a  natural  choice  is  linear  interpolation.  A  continuous  //-piecewise 
linear  function  is  also  a  continuous  /i-piecewise  linear  function.  That  is,  the  "coarse" 
finite  element  space  is  a  subspace  of  the  "fine"  finite  element  space  (see  [Ref.    1]). 
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llii:-  inc. ins  that  lineai-  interpolation  corresponds  to  the  "natural  embedding"  from 
the  coarse  finite  element  space  to  the  fine  element  space.  Its  advantages  are  its 
simplicity  and  its  linearity.  The  stencil  is  given  by 


i) 


B.      AMPLIFICATION  MATRIX 

Now  that  we  have  our  fine  grid  operators,  coarse  grid  operators,  and  intergrid 
transfer  operators,  we  conduct  analysis  using  the  DFT  (see  Appendix  A)  similar  to 
that  in  Chapter  IV  to  obtain  some  indication  as  to  how  a  two-level  multigrid  scheme 
using  these  operators  will  perform.  We  start  by  noting  that  the  process  u(  relaxing 
Dii  the  fine  grid,  obtaining  a  coarse  grid  correction,  and  relaxing  again  on  the  tine 
grid  can  be  represented  by  the  two-level  error  propagation  matrix,  PU2{(  '(i)I)l'x .  Note 
that,  while  h  as  a  superscript  identifies  grid  spacing,  U\  and  v<i  are  exponents.  Thus, 
P"1  and  P"2  represent  relaxation  U\  and  v2  times  on  the  fine  grid,  and  CG  is  the 
coarse  grid  correction.  Recall  that  coarse  grid  correction  is  the  following  series  of 
steps: 

Coarse  grid  correction  (two-grid  scheme) 


•  Relax  on  LhTh  =  /    on  iVl  to  obtain  an  approximation  v  . 

•  Compute  the  residual  rh  =  fh  —  Lhvh . 


Restrict  the  residual  rH  =  lj?rh. 


•   Solve  the  residual  equation  LH e     —  rH  on  0     to  obtain  an  approximation  to 
the  error  e    . 


Interpolate  the  error  el  =  I" e 


h   _    iH  CH 


•  Correct  the  approximation  obtained  on  iVl  with  the  error  estimate  obtained 
on  nH:  vh  <-  vh  +  e/l. 


rims.  CG  can  be  written  as 

CG=  I  -  lh{LH)-ltfL'\ 

where  /  is  the  identity  operator,  if?  and  I  fa  are  the  restriction  and  interpolation 
operators,  Lh  is  the  fine  grid  operator,  and  (LH)~l  is  the  coarse  grid  smoother.  Each 
of  the  operators  in  the  two-level  error  propagation  matrix  is  expanded  in  a  DFT 
(see  Appendix  A  and  Chapter  IV),  making  the  resulting  expansions  functions  of 
/.  ???  =  — y  +  1  :  y,  the  indices  of  the  expansion  in  the  frequency  domain.  These 
expansions  are  used  to  construct  an  amplification  matrix,  A(/,m),  which  is  also  a 
function  of  /  and  m.  The  largest  spectral  radius  of  this  matrix  over  the  values  that 
/.  //?.  have  on  the  coarse  grid,  i.e.,  l,m  —  —  ^  +  1  :  -j  (this  relationship  is  discussed 
later),  will  give  us  a  two-level  asymptotic  convergence  factor,  A.  This  factor, 
much  like  the  values  determined  for  the  relaxation  schemes,  will  give  a  bound  on 
how  well  the  two-level  scheme  can  be  expected  to  perform.  In  order  to  guarantee 
convergence,  A  <  1  is  necessary;  the  smaller  A  is,  the  better. 

In  Chapter  IV,  we  used  DFT  expansions  of  the  error  equations  derived  from  the 
relaxation  operators  to  determine  the  frequency  domain  coefficients  for  the  operators. 
We  now  use  the  same  technique,  again  using  error  equations,  to  determine  the  matrix 
entries  for  the  operators  in  the  two-level  scheme.  The  amplification  matrix  is  formed 
by  multiplying  together  these  operator  matrices, 

A(/,  m)  =  (PY2 (I  -  lkH{LH)-li1l!tk)(PhP  ,  (VI.3) 

where 

•  hh  is  the  4x4  matrix  for  the  fine  grid  discretization  operator. 

•  Ph  is  the  4x4  matrix  for  the  fine  grid  relaxation  operator. 

•  (Lw)_1  is  the  lxl  matrix  for  the  coarse  grid  relaxation  operator. 

•  1^  is  the  1x4  matrix  for  the  restriction  operator. 


• 


lit  is  the  4x1  matrix  for  the  interpolation  operator. 
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•   I  is  the  -1x1  identity  matrix,  (see  [Ref.  2]) 

We  go  through  the  details  oi  computing  the  matrix  entries  for  the  tine  grid  operator 
/.  ' :  results  fur  the  remaining  operators  are  then  presented. 

Prom  the  FVE  stencils  for  L  ,  we  have  the  operator  equation 


Ji(new)    _       ,  /i     li(old) 


2k,2j 


'2k, 2  j 


or 


li(new) 

V2A;,2./ 


<7384 


+(6«-ll)»Jg?-i 
I   +(32j  -  5)«J£  *,-_, 


(64i-5)^;2j      +(32j  +  5)^i+1 
+448jvjgf  +(64j  +  llK^f+1 

f  (64j  +  5)^t%3 


OK/i 


2 


/ 


V 


/l(0l(i) 


4-7U2A:+l,2j 

h(old) 


h(old) 


+(-i  +  f;>23LI    -i6iv23'  +(i  +  4i)«3s;Vi 


\ 


/ 


where  the  superscript  /i  indicates  a  vector  component  on  the  fine  grid,  (old)  and 
[new)  indicate  before  and  after  the  operator  acts,  and  subscripts  2A;  and  2j  are  the 
respective  fine  grid  row  and  column  indices.  These  indicial  values  are  used  because 
we  shortly  will  need  to  discuss  the  correspondence  between  the  fine  and  coarse'  grid 
points;  only  those  fine  grid  points  with  indices  2A;  and  2j  correspond  to  coarse  grid 
points  (with  indices  j  and  k).  Ordinarily  (see  [Ref.  2]),  the  consideration  of  indicial 
values  is  not  a  concern.  In  this  case,  however,  since  the  stencils  are  grid  location 
dependent,  and  since  the  analysis  of  the  amplification  matrix  must  apply  to  both 
the  fine  and  coarse  grids,  the  j  that  is  used  is  the  column  index  on  the  coarse  grid, 
requiring  that  2j  be  the  column  index  on  the  fine  grid.  Thus,  expanding  both  sides 
in  a  DFT,  we  have 
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E      E    ^c(i,m)  =     E      E     si((W-5)#(i,'«)«i 


/  =  --f  +  l  m  =  -f  +  1 


i=-£  +  im=_4  +  i 


/r' 


2-n-(l  +  m) 
JV 


+  (32j+5)V^C(/,m)e 
+(64j-llW£)C(/,m)e^ 
+448jV$C(l,m)  +  (64j  +  ll)V$C(l,m)e^ 

+(32j-5)l&)C(Z,m)e — 

+  (64j  +5)VjSC(Z,m)c-^) 


■  2ttI 


1277T7I 


+(-l+4i)V/JJC7(/,m)c 
-16jV;SC(/,m)  +  (H-4i)VlSC(/,m)Ci 
+  4i^C(/,m)e^ 


27rm 
N 


-i2tt/ 


where  2  <  2j  <  /V  —  2  on  a  grid  labeled  0  :  yV,  since  we  are  in  the  interior,  /,m 

-,  and 


2    ~  2 


0(1,111)  =  e   "   e    N 


When  we  transfer  information  from  the  fine  grid  to  the  coarse  grid,  we  must 
remember  that  not  all  of  the  freqencies  on  the  fine  grid  can  be  represented  on  the 
coarse  grid.  The  highest  frequency  that  can  be  resolved  on  any  grid  is  the  Nyquist 
frequency  (see  [Ref.  12]),  /  =  ^x-?  where  Ax  is  the  grid  spacing.  Hence,  on  the 
coarse  grid,  QH ,  Ax  —  -^  and  /  =  -j.  Therefore,  in  order  to  make  this  analysis 
germane  to  both  fine  and  coarse  grids,  we  must  rewrite  the  sums  in  the  expansions 
so  that  |/|,  |m|  <  -j.  This  can  be  done  as  follows: 


J2         £      Vi,mC(l,m)=     £         E      N,mC(l,m) 


(VI.4) 


l=-%  +1  m= 


+1 


:=-f+im=-f+i 


N  N  N  N  I 

+Vl+KmC(l  +  —,m)  +  Vl>m+KC(l,  m+-)  +  Vl+Njm+KC(l  +  y,m  +  -) 
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I'lms.  i  cu'iit  ing  the  DFT  expansion  of  the  operator  equation  to  consider  t  hose 
frequencies  that  can  be  represented  on  the  coarse  grid  and.  making  use  of  the  or- 
thogonality property  o\  the  complex  exponential  (see  Appendix  A),  we  can  equate 

individual  terms  of  the  sum,  and  t  hen  divide  i)\r  C(l,m)  f-  0  to  give 

k3 


vr 


£7384 


.     .     Af,\  '2Tl  .     ,  |  , .  >  l27T(t  +  T7l) 

64j    -5)VjSe-Tr+(32j  +  5)V;Sc- 


•  i  - 


i  a 


(o) 


+(64i-  u)^—  +448yv^  +  (64j  +  ii)\;:;;v 


+(32i-5)V/Sc" 


!2ir((  +  rn| 


.(.,) 


+  (64i  +  5)0 


O  A,' , 


■M.H 


T-(W,>Jlr  +  (-l+4i)^1 


(«)     -J 


•I") 


•H,^ 


(o) 


-16j^  +  (l+4j)^c^r+4j^c 


or 


,  factoring  out  Vj  °t, 


64j  -  5)e  "    +  (32j  +  5)e      -v 


flf384 


+(64  j  -  ll)e^~  +  448/  +  (64j  +  ll)e 


+(32 j  --  5)e  +  (64j  +  5)e     " 


l2  JTTTI 

N 


(4je  w    +  (-1  +  4j)e      -v 

%        i27rm  -i2n-<  \ 

16J +  (1 +4j)e   N     +4je^v-J 


r, 


(0) 


l, m  1 


and 


h 


(?384 
+(64;--  H)e 


(64j  -  5)eV   -  (32j  +  5)e 


i2jr(i  +  m) 


+  448j  +  (64j  +  ll)e" 


t2ir(l  +  m) 

■(32 j  -  5)e "       -  (64 j  +  5)e" 


i2tI 


CVK/ 


— 4je  w    +  (-1  +4j)e' 


l2TTTn 


16j  +  (l+4j)e—   -4je^-)]^\m, 
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V,         ,N 


g3S4 


(64j 


i27ri  i27r(f+™) 

5)e  w    -  (:y2j  +  5)e      " 


-(64 j  -  1  l)e-N-  +  448  j  -  (32  j  +  ll)e 


,  ■27r(i  +  m)  .271-1 
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t2TTTTt 
N 


atih  :2nl 

—  (4je  n    -  (-1  +  4j)e 
-16j-(l  +  2j)e—  +  4je 


_  i2irm 

N 

-i2irl 

N 


y(o) 


and 


since 


VK  ' 


.27r((+4) 


!r 


g384 
-(64j-ll)e 


64j  -  5)eV  +  (32j  +  5)e 


■  27T(i+m) 
N 


+  448j  -  (64j  +  ll)e" 


■  2TT(t+m)  I27TJ 

+(32j  -  5)e        "       -  (64j  +  5)e~~ 
^-(-4je—  -(-l+4j) 


2 
16j-(l+4j)e 


i27rm 


W  = 


i2si   i 
e  n  e 


■  2tt1 

—  e  ^ 


and 


4j/)e      n 

— t2ir(  \  1 

Aje    n    \ 

y(o) 

.27r(m+^) 

t27rm 

=  —  e  ^    . 

The  matrix  of  coefficients  for  the  fine  grid  operator  is  a  4x4  diagonal  matrix, 
Lk  =  (Aj,;)  (see  [Ref.  2]).  To  see  why  the  matrix  is  4x4,  recall  that  the  sum  in  the 
original  DFT  expansion  of  the  operator  equation  was  rewritten  so  that  the  indices 
satisfy  |/|,  |m|  <  — .  One  result  of  this  is  to  rewrite  each  of  the  individual  terms  of 
the  sum  as  a  combination  four  terms  (Equation  VI. 4).  That  is,  for  each  dimension 
of  the  problem,  represented  by  I  and  m,  there  are  two  components  of  each  individual 


term  in  the  sum,  which  are  indexed  by  |/|,  \m\  <  ^  and  -j  <  |/|,  \m\  <  y.  Thus,  the 
matrix  of  coefficients  is  a  22x22  matrix,  where  each  row  in  the  matrix,  say  the  z'th  row, 
comprises  the  coefficients  used  in  determining  the  corresponding  zth  component  in  the 
new  vector  as  a  linear  combination  of  components  of  the  old  vector.  In  this  case  the 
matrix  is  diagonal,  since  for  each  of  the  four  components  of  the  new  vector,  the  only 
contribution  from  the  old  vector  comes  from  the  corresponding  component.  In  other 
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ords,  the  only  contribution  to  V.^i      comes  from  V,%     .  So,  V/(n)  =  \jtl\'{").  where 


l+%,m 


l+f  ,7 


\ ''"'  .Hid  \  '"'  are  the  4x1  vectors  whose  components  have  indices  (/,  m),  (/  +  -r,  m), 


\ 


N 


I .  in  +  — ),  and  (/  +  —■ '"  +  -j).  The  I  (Mir  diagonal  entries,  A,-,,-,  are: 


A, 


^384 
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</384 
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QfAC/i     /        .      i2irl  i27rm 

2~  (4?e  "     _(_i+4j)e      w 
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</384 
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16;-(l+4;>— -4;V 


A  similar  process  must  be  applied  to  determine  each  of  the  other  matrices 
that  comprise  A(/,  m).  In  order  to  compute  the  entries  for  the  matrix  P/l,  recall  from 
Chapter  IV,  Equation  IV. 25,  that  we  can  write  the  error  equation,  enew  =  Pe°ld,  as 


—     (new)  OiKft    f     (iiew) 
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(new) 


{  —  /iv2k-l,2j-l  ~  av2k-\,2j  ~  L'V2k;2j-\ 


C/384 

njU2k,2j  +  l 


p  (old)  nJ°ld)         \ 


2k+l,2j         K-TV2k+\,2j  + 

\—nv2k_l2j  —  iv2k2j_1 


(old) 


(old) 
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where  A  =  32;  -5,5  =  64;  +  5,  C  =  64;  -  11,  D  =  448;',  F  =  64;  +  11,  F  = 
64;  -  5,  C  =  32;  +  5,  and  H  =  4;,  /  =  -1  +  4;,  J  =  - 16;',  K  =  1  +  4;',  L  =  4;.  Using 
similar  calculations  to  those  above,  we  have  the  four  entries,  n^;,  of  the  diagonal 
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In  similar  fashion,  the  entry  for  the  ixl  matrix  (Lw)    l,  corresponding  to  the 
coarse  grid  relaxation  operator,  is  given  by 


i2nll  +  m) 


]_oizk[_KHe^_LHe^ 


</3S4l 


i2n(t+m) 


+  ^r"  +  C^e^  +  £>"]  -  zf[H»e=W-  +  /"e^  +  JH] 


■here  ,t"  =  1 6j  -  5,  £"  =  32j  +  5,  C"  -  32  j  -11,  D"  =  224  j,  E"  =  32 j  +11.  /•,// 


32j  -  5,  G,/y  =  16j  +  5,  and  tf  w  =  2j,  /H  =  -1  +  2j,  ./"  =  -87,  A'w  =  1  +  2j,  Z."  =  2j. 
The  equation  for  the  conservation  restriction  operator,  Ih  ,  is 
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Let  Q  =  ^'-n  and  fi  =  \2 .  |,  and  expand  this  relation  in  a  DFT.  Considering 
those  frequencies  that  can  be  represented  on  the  coarse  grid  and,  making  use  of  the 
orthogonality  property,  we  can  equate  individual  terms  of  the  sum.  We  divide  by 
C(/,m)  =  e— ~e  n  ,  the  coarse  grid  analog  to  C(/,7n),  to  give  the  components  of 
the  matrix  of  coefficients  corresponding  to  the  restriction  operator,  ij^  =  (Ai)t).  This 
matrix  is  1x4  since  the  output  of  the  operator,  a  single  component  on  the  coarse  grid. 
is  a  linear  combination  of  22  components  on  the  fine  grid.  The  entries  in  the  matrix 
are  as  follows: 


Au  = 

Ai,a  = 

A,,3  - 

A!, 4  = 


1    +  CCS(-rr)     r  i2nm  „     %2nm 

N-  \Qe-—  +2  +  Re— 


2 


•1 


1 2  tt  tti  %2nm 

Qe     ~  +  2  +  Re~N~ 


^(-r)r.Qe-*p+2_Re 


\2-nm 

N 


and 


i-cosm) 


i2irm  t2-rrm 

—Qe     w    4-  2  —  /te   w 


SI 


The  equation  for  the  injection  restriction  operator  is  vjfj  —  t>2fc,2i  which,  using 
similar  calculation,  results  in 

V"n   =    Km  +  Vl+f,m  +   Kn+i   +   Vl+f  ,m+f  ' 
SO  that   l£   =   [1        1        1        1]. 

The  linear  interpolation  operator  gives  rise  to  the  following  relations: 

h  _       H 

V2k,2j      ~       Vk,ji 


V 


'2k+l,2j      —  r, 

2fc,2j  +  l  r> 


■k,j    T    "fc+l,J  +  l 

';2fc+l,2j  +  l       —  9 

Expanding  these  in  a  DFT,  considering  the  coarse  grid  frequencies,  and  making  use 
of  the  orthogonality  property,  we  can  equate  terms  of  the  sums  and  divide  by  C(l,m) 
to  give 

VlHm      =       Km  +  Vl+f,m  +  Km+f   +  Vl+%  ,m+f ' 

l27Ti  -1271-i 

yU    (e    N      +   e      N       )         _        T/Zi  T//l  ,      T//l  T/A 

.      i2TTTn  —  i2nm  . 

H  (e   n     +  e    at     )  .  .  .  , 

V/,m g ~  ''7n  +   Vf  ,m  +   V+f   +   Vf  ,m+f  '         and 

i27rl  +  m  -i27r(l  +  m) 

wH    \e       N         +  e  "  )  I/i        ,     T//1  ,     vh  ,     T/A 

l//,m 2  -       V/,m  +   Vf  ,m  +    Vm+f   +   ^+f  ,m+f " 

The  components,  0t  j ,  of  the  matrix  of  coefficients,  1^,  are  found  by  solving  the  above 
system  of  four  equations  for  the  fine  grid  components,  such  that  Vh  —  \hHVH .  The 
matrix  1^  is  4x1  since  the  output  of  the  operator,  the  22  components  of  the  fine  grid 
vector,  is  the  result  of  interpolation  operator  acting  on  the  single  component  of  the 
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coarse  grid  vector.  That  is, 


"  elf] ' 

e2ll 

03,1 

1 

1 

.   a>>    . 

1  +  cos^f)  +  cos(^)  +  cos(!^) 

1  -  cos(2-^)cos(^)  -  cos(^^) 
1  +  cos(^)  -  cos(2-Jr)  +  cos(2-^) 

1  -  cos(^)  -  cos(^)  +  co5(^±^) 

Thus,  we  have  all  the  elements  to  construct  the  amplification  matrix,  A(/.  m  ). 
i  Equation  VI. 3).  We  now  have  the  tools  to  analyze  the  expected  perfomance  of  a 
two-level  scheme. 

C.       NUMERICAL  RESULTS 

We  now  conduct  a  number  of  numerical  experiments  to  investigate  the  perfor- 
mance of  the  multigrid  method,  again  using  Matlab.  First  we  examine  the  behavior 
of  the  amplification  matrix,  A(/,?/i),  and  the  largest  spectral  radius  of  the  family  of 
matrices  A(/,ra)  for  all  /,m.  We  then  present  the  results  of  the  repeated  use  of  the 
V-cycle  to  solve  Equation  IV. 2,  as  compared  with  the  iterative  solutions. 

In  order  to  make  use  of  the  amplification  matrix,  we  must  choose  a  relaxation 
method,  and  the  number  of  times  to  relax  on  the  way  down  and  on  the  way  back  up, 
i.e..  we  must  choose  u\  and  u-i  in  Equation  VI. 3.  We  experiment  first  with  Gauss- 
Seidel  relaxation  and  three  different  combinations  of  [v\^v-i)  V-cycles:  a  (0,  1)  cycle, 
a  (2,1)  cycle,  and  a  (10,10)  cycle  (recall  that  a  =  1,  meaning  fully  implicit  time 
stepping,  and  the  time  step  g  =  h).  The  first  two  types  of  cycles  are  fairly  common 
(see  [Ref.  3]),  while  the  latter  is  included  for  reference.  The  information  presented  in 
Tables  VIII,  IX,  and  X  indicates  the  asymptotic  convergence  factor,  A,  for  specific  grid 
sizes,  grid  positions,  and  type  of  restriction  operator,  either  conservation  or  injection. 
Linear  interpolation  is  used  as  the  interpolation  operator  in  all  of  these  experiments. 
Tables  XI  and  XII  indicate  the  asymptotic  convergence  factors  for  horizontal  and 
vertical  line  relaxation  and  a  (2, 1)  cycle. 

Table  VIII  indicates  the  asymptotic  convergence  factors  that  result  from  a 
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(0,1)  cycle,  Gauss-Seidel 

Conservation 

Injection 

2j 

Ar  =  16 

N  =  32 

N  =  64 

.V  =  16 

N  =  32 

N  =  64 

2 

.9174 

.9651 

.9843 

.9397 

.9710 

.9858 

2 

.8308 

.9111 

.9544 

.9144 

.9553 

.9771 

/V  -  2 

.7472 

.8595 

.9259 

.8929 

.9422 

.9700 

fable  \  III.  Asymptotic  Convergence  Factors  for  the  (0,  1)  Cycle  Using  Gauss-Seidel 
Relaxation. 

single  iteration  sweep  on  the  fine  grid  coupled  with  the  course  grid  correction  for 
both  the  conservation  and  injection  restriction  operators.  In  general,  the  conservation 
restriction  operator  seems  to  indicate  better  performance,  since  its  convergence  factors 
are  lower  than  those  for  the  injection  operator  (recall  that  the  asymptotic  convergence 
factor  is  the  largest  spectral  radius  of  the  family  of  amplification  matrices,  which 
comes  from  the  DFT  expansion  of  each  component  in  the  two-level  error  propagation 
matrix).  Similar  to  the  analysis  on  relaxation  schemes  in  Chapter  IV,  the  convergence 
factors  decrease  with  radial  distance,  indicating  an  improvement  in  performance,  but 
increase  with  grid  size.  Even  though  the  convergence  factors  for  all  three  grid  sizes  are 
less  than  one,  indicating  that  convergence  is  very  likely  (but  not  necessarily  certain, 
since  these  factors  tell  us  nothing  about  what  is  going  on  at  the  boundaries),  at  least 
some  of  the  factors  for  all  three  grids  are  relatively  high,  indicating  that  convergence 
may  be  quite  slow  with  the  (0,1)  cycle.  A  point  of  interest  is  to  compare  these 
results  with  those  in  Table  IV,  which  indicate  the  maximum  ratio  of  new  to  old  error 
components  for  Gauss-Seidel  relaxation.  One  might  expect  the  results  in  Table  VIII 
to  be  universally  better  than  those  in  Table  IV,  however  this  is  not  the  case.  In 
fact,  on  a  grid  of  size  A^  =  32,  the  largest  asymptotic  convergence  factors  for  both 
the  conservation  and  the  injection  restriction  operators  are  larger  than  the  largest 
ratio  for  just  a  single  Gauss-Seidel  sweep  on  the  fine  grid.  In  other  words,  a  single 
Gauss-Seidel  relaxation  sweep  is  predicted  to  do  better  than  a  (0, 1)  cycle.  This  is  a 
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(2,1)  cycle,  Gauss-Seidel 

Conservation 

Injection 

2j 

N  =  16 

N  =  32 

N  =  64 

n  =  16 

N  =  32 

N  =  64 

2 

.8223 

.9133 

.9574 

.8423 

.9189 

.9589 

N 
2 

.7375 

.8572 

.9254 

.M  16 

.8987 

.9474 

N  -2 

.6623 

.8083 

.8977 

.7914 

.8861 

.9404 

liable  IX.  Asymptotic  Convergence  Factors  for  the  (2,1)  Cycle  Using  Gauss-Seidel 
Relaxation. 

facl  worth  remembering  when  the  multilevel  solution  is  analyzed. 

Table  IX  indicates  the  asymptotic  convergence  factors  for  a  (2,1)  cycle;  as 
expected,  the  factors  are  lower  than  the  corresonding  factors  for  the  (0, 1)  cycle,  and 
are  universally  lower  than  the  ratios  associated  with  a  single  Gauss-Seidel  sweep  as 
indicated  in  Table  IV.  However,  another  comparison  is  to  raise  the  ratios  in  Table  IV 
to  the  third  power,  corresponding  to  three  relaxation  sweeps  on  the  fine  grid.  In  other 
words,  since  for  a  (2, 1)  cycle  there  are  two  relaxation  sweeps  on  the  fine  grid  before 
the  course  grid  correction,  and  one  afterward,  for  a  total  of  three  relaxation  sweeps 
on  the  fine  grid,  it  seems  reasonable  to  compare  the  performance  of  a  (2, 1)  cycle  to 
three  iterations  of  the  Gauss-Seidel  method,  represented  by  raising  the  Gauss-Seidel 
factor  to  the  third  power.  For  example,  in  Table  IV  for  N  =  32,  j  =  1,  the  ratio 
is  .9567  which,  when  raised  to  the  third  power  is  .8756.  This  value  is  lower  than 
either  of  the  values  in  Table  IX  for  N  =  32,  2j  =  2.  This  also  holds  true  for  the 
asymptotic  convergence  factor  associated  with  N  =  32  and  2j  =  ~  (this  phenomenon 
does  not  occur  for  N  =  32,  2j  =  N  —  2,  nor  for  tV  =  16).  Thus,  the  overall  predicted 
performance  of  three  Gauss-Seidel  relaxation  sweeps  for  N  =  32  is  better  than  the 
predicted  performance  of  a  Gauss-Seidel  (2,  1)  cycle. 

Similar  to  the  information  in  Table  VIII,  the  performance  predictions  in  Ta- 
ble IX  improve  with  radial  distance,  and  worsen  with  increase  in  grid  size.   Also,  as 
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(10,10)  cycle,  Gauss-Seidel 

Conservation 

Injection 

2j 

/V  =  16 

N  =  32     N  =  64 

N  =  16 

N  =  32 

N  =  64 

2 

.3243 

.5716         .7567 

.3322 

.5751 

.7579 

N 
2 

.2678 

.5104         .7118 

.2948 

.5351 

.7288 

N  -  2 

.2374 

.4795         .6898 

.2837 

.5256 

.7226 

Table  X.  Asymptotic  Convergence  Factors  for  the  (10,  10)  Cycle  Using  Gauss-Seidel 
[Relaxation. 

above,  the  conservation  restriction  appears  to  predict  better  performance  than  injec- 
tion. The  factors  are  still  relatively  high,  especially  for  the  N  =  64  grid  size,  so  that 
convergence  will  likely  be  slow. 

The  (10,10)  cycle  is  included  for  reference  so  that  the  performance  of  the 
two-level  cycle  can  analyzed  using  a  large  number  of  iteration  sweeps  on  the  fine  grid. 
The  characteristics  of  the  (10,  10)  cycle  are  similiar  to  those  of  the  (0,  1)  cycle  and  the 
(2.  1 )  cycle;  the  performance  improves  with  radial  distance  and  worsens  with  increase 
in  grid  size.  The  asymptotic  convergence  factors  are  smaller  than  for  previous  cycles, 
but  this  is  expected  since  a  good  deal  more  work  is  represented  by  the  (10,  10)  cycle. 
However,  using  a  comparison  similar  to  that  in  the  discussion  of  Table  IX,  in  Table  IV 
for  .'V  =:  32,  j '•  =  1,  the  ratio  is  .9567  which,  when  raised  to  the  20th  power  is  .4126. 
Here,  as  in  the  case  of  a  (2,1)  cycle,  this  value  is  lower  than  either  of  the  values 
in  Table  X  for  N  =  32,  j  =  1.  This  holds  true  for  all  the  factors  associated  with 
.V  =  32,  but  is  not  true  for  yV  =  16.  Thus,  as  before,  the  predicted  performance 
of  Gauss-Seidel  relaxation  for  TV  =  32  is  better  than  the  predicted  performance  of  a 
Gauss-Seidel  two-level  cycle. 

Tables  XI  and  XII  indicate  the  asymptotic  convergence  factors  for  horizontal 
and  vertical  line  relaxation  (2, 1)  cycles.  Similar  to  the  results  of  the  analysis  on  the 
horizontal  and  line  relaxation  schemes  in  Chapter  IV,  both  line  relaxation  schemes  for 
the  (2,  1)  cycle  indicate  better  performance  than  that  for  the  Gauss-Seidel  method; 
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(2,1)  cycle.  Horizontal  Line  Relaxation 

Conservation 

Injection 

2j 

N=  =  16 

N  -=  32 

.V  :=  64 

.V  =  16 

N  ■■=  32 

N  =  64 

2 

.6784 

.8262 

.9097 

.6938 

.8310 

.9111 

N 
2 

.6210 

.7836 

.8839 

.6793 

.8203 

.9046 

A/ -2 

.5629 

.7408 

.8580 

.6647 

.8096 

.8981 

Table  XI.  Asymptotic  Convergence  Factors  for  the  (2,  I  )  Cycle  Using  Horizontal  Line 

Relaxation. 


(2,1)  cycle,  Vertical 

Conservation 

Injection 

2j 

N  =   16 

N  =  32 

.V  =  64 

N  =  16 

N  =  32 

A   64 

2 

.7258 

.8562 

.9265 

.7429 

.8613 

.9279 

N 
2 

.6326 

.7876 

.8851 

.6928 

.8247 

.9059 

N-2 

.5685 

.7  127    .8586 

.6723 

.8119 

.8988 

Table  XII.  Asymptotic  Convergence  Factors  for  the  (2,1)  Cycle  Using  Vertical  Line 
Relaxation. 

horizontal  line  relaxation  is  the  better  of  the  two.  One  difference  is  that  in  Tables  V 
and  VI,  the  predictions  for  horizontal  line  relaxation  worsen,  and  the  predictions  for 
vertical  line  relaxation  improve,  with  radial  distance,  whereas  in  Tables  XI  and  XII 
the  predictions  for  both  line  relaxation  schemes  improve  with  radial  distance.  Using 
calculations  similar  to  those  in  the  discussions  of  Tables  IX  and  X,  the  values  in 
Tables  V  and  VI  are  raised  to  the  third  power  and  compared  with  the  factors  in 
Tables  XI  and  XII.  The  result  is  that  the  predicted  performance  of  three  sweeps  of 
the  horizontal  line  relaxation  is  better  than  that  for  the  (2, 1)  cycle  for  both  N  =  16 
and  N  —  32.  The  predicted  performance  of  three  sweeps  of  the  vertical  line  relaxation 
is  better  than  that  for  the  (2, 1)  cycle  only  for  N  =  32.  These  results  are  consistent 
witli  those  of  the  Gauss-Seidel  relaxation  lor  all  of  the  (1/1,1/2)  cycles  tested.  That 
is,  for  A/  =  32,  the  predicted  performance  of  ordinary  iteration  is  better  than  the 
predictions  for  any  of  the  corresponding  two-level  schemes  considered  here. 
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In  the  initial  stages  of  using  V-cycles  to  solve  Equation  IV. 2  at  a  single  time 
step,  we  experiment  with  the  use  of  conservation  restriction  and  linear  interpolation. 
The  result  is  that  this  combination  does  not  result  in  a  coarse  grid  correction.  That 
is.  the  effect  of  applying  the  two-level  scheme  is  to  generate  an  approximation  that 
is  not  as  good  as  that  obtained  by  relaxation  on  the  fine  grid  alone.  Therefore,  we 
experiment  with  the  combination  of  the  injection  restriction  operator  and  the  linear 
interpolation  operator,  with  the  result  that  this  combination  does  produce  a  coarse 
grid  correction.  A  correction  scheme  using  V-cycles  with  these  intergrid  transfer 
operators  is  implemented  to  solve  at  a  single  time  step,  and  then  the  solution  is 
stepped  in  time  until  a  "steady"  solution  was  obtained. 

Multigrid  V-eycle  solutions  to  Equation  IV. 2  are  computed  on  grids  /V  =  8,  16, 
and  32,  using  (2.  1)  cycles,  Gauss-Seidel  relaxation,  a  —  1  (fully  implicit  time  step- 
ping), and  time  step  g  =  /?,,  where  the  initial  temperature  distribution  and  boundary 
conditions  are  the  same  as  in  the  iterative  solution.  We  again  use  the  discrete  energy 
norm  (Equation  V.l)  to  measure  the  accuracy  of  these  solutions, 

HIZ^HI  =  {LhDh,Dh)*, 

where  {•,■)  denotes  the  Euclidean  inner  product,  Lh  is  the  operator  defined  in  Equa- 
tion IV. 2,  and  Dh  —  Th  —  T*  approximates  the  discretization  error,  where  Th  is  the 
solution  to  Equation  IV. 2  on  grid  h,  and  T*  is  the  "exact"  solution  (i.e.,  solution  of 
Equation  IV. 2  with  /V  =  64  sampled  on  grid  h)  (see  [Ref.  11]).  The  results  presented 
in  Table  XIII  are  almost  identical  to  those  obtained  for  the  iterative  solution  (see 
Chapter  V). 

As  another  measure  of  how  well  the  iterative  solutions  match  the  multilevel 
solutions,  the  discrete  L2  norm  of  the  solution  differences,  d'\ 

\\dk\\   =   ±(dh,dh)L2, 

is  computed  for  each  of  the  above  grid  sizes  and  time  value.  The  results  are  presented 
in  Table  XIV;  while  the  norm  of  the  differences  is  generally  small,  it  increases  with 
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i  N  =  8 

V  =  16 

N  -=  32 

Initial  1  ime  step 

43.07 

23.97 

6.87 

Time  =  "\  second" 

42.99 

23.98 

6.87 

Steady  state                  42.99        23.98 

6.87 

; 

fable  XI 11.  Discretization  Errors  for  the  Multilevel  Solution. 


.'V  =  8 

N  ----  16 

.V  ==  32 

Initial  time  step 

s.diiE  It, 

1.74E-15 

3.95E-15 

Time  =  "1  second'1 

5.68E-16 

1.53E-15 

4.64E-15 

Steady  state 

3.02E-15 

1.02E-14 

5.66E-14 

fable  XIV.  Norms  of  Differences  Between  Iterative  and  Multilevel  Solutions. 

grid  size.  This  might  suggest  a  normalization  problem,  however,  the  factor  of  -^  in  the 
discrete  L2  norm  is  used  specifically  to  avoid  this  type  of  problem.  More  investigation 
is  required  to  determine  why  there  is  not  better  agreement  between  the  iterative  and 
multilevel  solutions. 

Figures  30  and  31  indicate  the  number  of  iterations  per  time  step  for  the  it- 
erative solutions  and  the  number  of  V-cycles  per  time  step  for  the  multigrid  solution 
for  N  —  32.  It  is  important  to  note  that,  while  both  processes  require  a  decreasing 
amount  of  work  for  each  time  step  out  to  about  150  steps,  the  multilevel  process  re- 
quires about  \  more  time  steps  than  the  iterative  process  to  reach  a  steady  solution: 
the  multilevel  process  requires  one  V-cycle  per  time  step  for  time  step  >  150.  In 
addition.  Figures  32  and  33  illustrate  how  t  lit-  norm  of  the  difference  between  suc- 
cessive solutions  behaves  as  the  solution  is  stepped  in  time.  As  would  be  hoped,  the 
difference  between  successive  solutions  decreases  with  time  stepping.  However,  as  is 
evident  from  these  latter  figures,  the  time  stepping  process  stalls  somewhat  before 
the  steady  solution  is  reached  (the  tolerance  used  is  machine  e  =  2.22E  —  16).  While 
the  number  of  time  steps  to  reach  a  steady  solution  on  grids  N  =  8, 16  (not  shown) 
is  about  the  same  for  both  the  iterative  and  multilevel  solutions,  on  grid  N  =  32, 
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as  already  indicated,  the  number  of  time  steps  needed  for  the  multilevel  process  to 
reach  steady  solution  is  about  ~  longer  than  that  required  for  the  iterative  solution. 
Even  so,  the  two  processes  on  all  three  grid  sizes  stall  at  about  the  same  time  step. 


Iterations  per  time  step 


200 


250 


Figure  30.  The  number  of  iterations  per  time  step  for  N  =  32;  170  time  steps  needed 
to  reach  a  steady  solution. 

We  conclude  our  analysis  of  the  multilevel  solution  with  a  discussion  of  how 
much  computational  work  must  be  done  to  reach  a  solution.  For  this  purpose,  we 
pick  a  time  step,  for  both  methods,  for  which  the  time  stepping  process  has  stalled. 
For  /V  =  8,  use  time  step  s  =  50;  for  N  =  16,  use  time  step  s  =  80;  and  for  N  =  32, 
use  a  time  step  of  s  =  150.  The  norm  of  the  difference  between  successive  solutions 
for  these  time  steps  is  on  the  order  of  10-15.  We  use  calculations  similar  to  those  in 
[Ref.  3]  in  computing  the  cost.  If  we  consider  a  work  unit,  WU,  to  be  the  cost  of 
one  relaxation  sweep  on  the  finest  grid,  then  we  can  estimate  how  many  work  units 
are  associated  with  each  V-cycle.  Since  we  are  using  a  (2,1)  cycle,  relaxation  occurs 
three  times  on  each  level;  one  relaxation  sweep  on  the  grid  Qph  requires  p~d  of  a  work 
unit,  where  d  is  the  dimension.  Adding  these  costs,  and  using  the  geometric  series  as 
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Figure  31.  The  number  of  V-cycles  per  time  step  for  N  =  32;  230  time  steps  needed 
to  reach  a  steady  solution. 

an  upper  bound,  we  have 


.->  —  nd 


V-cycle  Cost  =  3{1  +  Td  +  ■■■  +  2_nrf}WU  < 


1  —  2~d 


WU. 


Thus,    for  the  two-dimensional  case  we  have  Cost   =   4WU.     A  comparison  of  the 
computational  cost  for  the  various  grid  sizes  is  given  in  Table  XV. 


N  =  S 

N  =  16 

N  =  32 

Time  step  (s) 

50 

80 

150 

"Elapsed"  time  (jj) 

6.25 

5.00 

4.69 

Number  of  V-cycles 

1246 

11)7(1 

13281 

Total  Cost  (WU) 
V-cycle  cost  (WU) 

4984 

16304 

53124 

Iteration  cost  (WU) 

6245 

25631 

103185 

Table  XV.  Computational  Cost  for  the  Iterative  and  Multigrid  Solutions  (Not  Includ- 
ing the  Cost  of  the  Intergrid  Transfers). 
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Norm  of  difference  between  successive  time  step  solutions 


250 


Figure  32.   Norm  of  difference  between  successive  iterative  solutions  for  N  =  32;  170 
time  steps  needed  to  reach  a  steady  solution. 

The  cost  of  implementing  the  intergrid  transfers  is  not  included  in  the  com- 
putations presented  in  Table  XV.  In  considering  the  amount  of  work  represented 
by  one  work  unit,  counting  "adds"  and  "multiplies"  each  as  one  floating  point  op- 
eration (flop),  the  number  of  flops  for  one  relaxation  sweep  on  a  grid  of  size  N  is 
approximately  4  IN'2.  Injection  restriction  requires  no  arithmetic;  linear  interpola- 
tion requires  approximately  ^-  flops.  Thus,  the  cost  of  the  intergrid  transfers  for 
the  two-level  scheme  is  about  (^-/AIN2)WU  =  ^WU .  Again  using  the  geometric 
series  as  an  upper  bound,  we  have  that  the  cost  of  the  intergrid  transfers  is  given  by 

Transfer  Cost  =  — -{1  +  2~d  +  •  •  •  +  2_nd}WU  <  -^-WU. 
.328  246 

Thus,  we  update  the  information  in  Table  XV  and  present  the  results  in  Table  XVI. 
While  not  dramatic,  the  use  of  multigrid  does  result  in  a  savings  of  compu- 
tational cost.   The  cost  savings  increase  with  grid  size;  the  savings  on  a  grid  of  size 
/V  =  8  is  only  20%,  on  a  grid  of  size  N  =  32  the  savings  is  almost  50%.    In  consid- 
ering these  results,  several  questions  need  to  be  answered.    For  example,  are  these 
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Norm  ot  difference  Detween  successive  time  step  solutions 


.  50 


Figure  33.  Norm  of  difference  between  successive  multilevel  solutions  for  A  =  32;  230 
time"  steps  needed  to  reach  a  steady  solution. 


iV  =  8 

N  =  16 

N  =  32 

Total  Cost  (WU) 

V-cyclecost  (WU) 

4992 

16354 

53286 

Iteration  cost  (WU) 

6245 

25631 

103185 

Table  XVI.  Computational  Cost  for  the  Iterative  and  Multigrid  Solutions,  Including 
i  he  ( 'ost  of  the  Integrid  Transfers. 

results  consistent  with  the  local  mode  analysis  and  the  information  that  comes  from 
the  amplification  matrices?  Does  the  multigrid  solution  show  "typical"  performance; 
i Iocs  the  process  work  as  one  would  hope? 

In  general,  the  results  are  consistent  with  the  local  mode  analysis,  which 
predicts  that  convergence  will  be  slow.  The  prediction  of  the  amplification  matri- 
ces is  difficult  to  apply,  since  the  asymptotic  convergence  factors  apply  to  two-level 
schemes  and  not  to  V-cycles.  Nevertheless,  the  convergence  factors  did  not  predict 
that  the  two-level  scheme  would  provide  improvement  over  the  iterative  process  on 
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grid  size  N  =  32,  and  therefore  these  results  are  somewhat  of  a  pleasant  surprise. 
In  other  words,  multigrid  performed  better  than  might  have  been  expected  based  on 
t  he  two-level  convergence  factors.  However,  typical  multigrid  performance  reaches 
convergence  (to  the  level  of  truncation  error,  see  [Ref.  3]  and  [Ref.  2])  in  just  a 
few  V-cycles.  Thus,  the  multilevel  solution  developed  here  is  disappointing,  since 
convergence  requires  a  very  large  number  of  V-cycles. 

There  are  no  obvious  causes  for  the  disappointing  performance  of  the  multi- 
level solution  of  Equation  II. 2.  One  of  the  problems  may  be  the  use  of  cylindrical 
coordinates,  which  results  in  a  radial  bias  in  the  FVE  stencils.  The  discretization  of 
i  lie  time  derivative  using  finite  differences,  while  using  the  FVE  method  discretize 
I  lie  spatial  derivatives,  may  also  have  some  impact  on  the  solution  process.  It  is  not 
clear,  however,  that  either  of  these  affects  the  multilevel  solution  any  differently  than 
it  affects  the  iterative  solution.  The  choice  of  relaxation  scheme  is  a  major  factor 
in  determining  the  convergence  of  the  process;  perhaps  the  use  of  line  relaxation,  or 
some  other  type  smoother,  will  improve  the  performance.  Again,  however,  it  is  not 
certain  that  this  will  result  in  an  improvement  of  the  multigrid  performance  relative 
to  the  iterative  process.  The  determination  of  the  coarse  grid  and  intergrid  transfer 
operators  is  also  something  of  a  concern.  This  is  one  area  that  clearly  affects  the 
performance  of  the  multilevel  solution  as  compared  to  the  iterative  solution.  The  fact 
that  the  predicted  performance  of  iteration  sweeps  alone  is  better  than  the  predicted 
performance  of  the  two-level  schemes  emphasizes  the  need  to  consider  the  use  of  dif- 
ferent coarse  grid  and  intergrid  transfer  operators.  Finally,  the  point-source  problem 
includes  a  singularity  at  the  boundary.  Most  of  the  analysis  has  been  done  for  interior 
points,  and  the  solution  results  do  not  drastically  differ  from  what  is  predicted  by  the 
analysis.  Neverthless,  the  treatment  of  the  boundaries,  particularly  for  this  problem, 
may  provide  a  means  to  improve  the  solution  process. 
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VII.  CONCLUSION 


I 'lie  purpose  ol  t  his  work  is  to  apply  the  FVE  met  I  km  1  to  discrel  ize  the  axisym- 
metric  heat  equation,  and  implement  multigrid  in  solving  it.  In  addition  to  using  the 
FVE  method  to  discretize  the  spatial  derivatives,  finite  differences  arc  applied  to  the 
time  derivatives,  resulting  in  a  time-stepping  scheme  that  makes  use  of  a  weighted 
average  oi  the  current  and  future  time  steps.  The  use  of  cylindrical  coordinates  results 
in  a  discretization  stencil  that  has  a  radial  bias,  the  effect  of  which  is  evident  in  the 
analysis  ol  various  relaxation  schemes.  Gauss-Seidel  is  the  relaxation  scheme  chosen 
[or  implementation  in  the  iterative  and  multilevel  solution  processes.  Due  to  numer- 
ical anomalies  encountered  in  computing  the  iterative  solution,  the  weighted-average 
time-stepping  scheme  becomes  a  fully  implicit  backward  time-stepping  scheme.  Its 
use  requires  more  computational  work  than  the  weighted-average  scheme,  but  its  use 
is  required  to  achieve  a  physically  meaningful  solution.  The  specifics  ol  the  multilevel 
technique  are  then  presented,  including  the  development  of  the  coarse  grid  and  inter- 
grid  transfer  operators,  and  the  predictions  of  expected  performance  for  the  two-level 
scheme.  The  multilevel  solution  is  then  implemented;  its  results  are  analyzed  and 
compared  to  the  results  of  the  iterative  solution  process.  The  results  of  our  work  are 
somewhat  disappointing  in  that  we  are  not  yet  able  to  achieve  the  hoped-for  level  of 
success.  There  are,  however,  a  number  of  positive  results  that  come  from  this  work, 
specifically,  the  application  of  the  FVE  method  to  a  problem  in  cylindrical  coordi- 
nates and  the  use  of  Maple  software  to  compute  the  necessary  integrals.  Additionally, 
we  know  that  the  usual  multilevel  techniques  do  not  produce  the  expected  results  and 
therefore  there  is  ample  room  for  further  study. 

There  are  several  possibilities  for  future  work.  One  specific  question  that  has 
not  been  addressed  is  the  determination  of  the  consistency  of  the  discretization  oper- 
ator.  Numerical  results  are  useful  in  attempting  to  verify  this  result,  but  consistency 
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should  be  proven  analytically.  In  conjunction  with  this  idea,  there  may  be  a  require- 
ment to  review  the  assumption  that  the  solution  can  be  adequately  represented  by 
a  collection  of  continuous,  piecewise  linear  basis  functions.  There  is  no  evidence  as 
yet  that  would  require  this  review,  but  it  is  possible.  Additionally,  we  are  assum- 
ing a  uniform  step  size  that  applies  to  both  the  r  and  z  directions.  One  possible 
modification  to  this  approach  is  to  implement  a  non-uniform  grid,  perhaps  one  that 
results  in  control  volumes  that  are  the  same  for  every  point  on  the  grid.  The  point  of 
this  modification  would  be  to  eliminate  the  radial  bias  in  the  discretization.  For  the 
discretization  of  the  time  derivatives,  the  requirement  to  produce  a  physically  mean- 
ingful solution  has  resulted  in  a  fully  implicit  time-stepping  scheme.  While  there  may 
l>e  another  method  to  discretize  the  time  derivatives,  it  will  have  to  meet  the  same 
criterion  and,  therefore,  may  not  be  much  of  an  improvement. 

All  of  the  above  considerations  apply  generally  to  the  problem,  and  not  specifi- 
cally to  multigrid.  There  are  at  least  three  areas  in  which  multigrid  performance  may 
be  improved:  intergrid  transfer  operators,  coarse  grid  operator,  and  treatment  of 
boundary  conditions.  The  intergrid  transfer  operators  used  are  the  simplest  available 
that  result  in  coarse  grid  correction  and,  consequently,  choosing  different  operators 
may  improve  performance.  Other  possibilities  include  the  use  of  a  restriction  operator 
that  is  based  on  the  piecewise  linearity  of  the  basis  functions,  or  perhaps  an  inter- 
polation operator  that  is  based  on  enforcing  conservation.  In  other  words,  choose 
intergrid  operators  that  satisfy  the  variational  condition,  ifa  =  c(lff),  for  c  some 
constant  (see  [Ref.  3]  and  [Ref.  1]).  The  choice  of  the  coarse  grid  operator  may 
need  to  be  made  in  conjunction  with  these  intergrid  operators.  That  is,  despite  the 
computational  complexities  that  would  be  introduced,  choose  a  coarse  grid  operator 
such  that  the  Galerkin  condition  Lh  =  I  fa  LH  iff  (see  [Ref.  3]  and  [Ref.  1])  is  sat- 
isfied. Additionally,  special  treatment  of  the  boundaries  may  be  required,  since  the 
boundary  of  the  point-source  problem  has  a  singularity  at  the  origin  (see  [Ref.  2]). 

Finally,  the  consideration  of  the  point-source  problem  is  the  first  step  in  ap- 
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plying  the  FVE  discretization  and  multilevel  solution  to  the  Navier-Stokes  equation 
i  hat  arises  from  the  molten-pool  problem.  'Tie  intricacies  oi  the  molten-pool  problem 
present  several  challenges.  Included  are  the  requirement  to  solve  a  system  ol  three 
PDEs  in  three  unknowns  and  the  requirement  to  keep  track  of  the  moving  phase- 
change  boundary  as  the  molten  pool  expands.  Additionally,  high  flow  velocities  and 
small  local  length  scales  result  when  convection  is  vigorous  ([Ref.  13]),  causing  fur- 
ther numerical  complications.  It  is  anticipated  that  the  Full  Approximation  Scheme 
(  FAS )  (see  [Ref.  2]  and  [Ref.  1  ] )  can  be  used  to  efFecl  ively  treat  the  non-linear  nature 
of  the  problem,  and  that  the  Fast  Adaptive  Composite  Grid  (FAC)  method  (sec5  [\\vl\ 
I J  will  be  useful  in  overcoming  difficulties  that  arise  in  the  geometry  of  the  problem. 
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APPENDIX  A.  THE  DISCRETE  FOURIER 

TRANSFORM 

Below  we  briefly  list  the  defiiiil  ion  and  some  of  I  he  propert  ies  of  t  he  Discrete 
Fourier  Transform.  DFT.  For  an  exhaustive  treatment  sec  [Ref.  10];  what  appears 
in  this  appendix  closely  follows  that  treatment. 

The  DFT  is,  in  at  least  one  interpretation,  a  way  to  discretely  approximate 
the  Fourier  transform.  For  example,  assume  we  are  working  on  a  temporal  or  spatial 
domain  [  —  7,  "f]  Wltn  grid  spacing  Ax  =  A/N  and  grid  points  x3  =  jAx,  and  its 
associated  frequency  domain  [—  — ,  -j]  with  grid  spacing  AlJ  =  2ir/A  and  grid  points 
.j.„  =  mAw,  where  j,  m  =  —^j  +  l  '■  t-  Suppose  r,  denotes  the  sampled  values,  ri.r.i. 


01  a  tui 


notion  defined  on  the  domain  [—7,  4],  and  that  v(u)m)  is  the  Fourier  transform 


of  v  at  the  frequency  grid  points,  ujm,  where  the  Fourier  transform  is  given  by 

/. , 
v{x)e~t2iruxdx. 
-00 

Using  the  Trapezoidal  Rule  to  approximate  the  integral,  we  have 


Vje     n 


where  m  =  — y  +  1  :  y.  Given  the  set  of  N  sample  values  of  Uj,  the  DFT  comprises 
the  N  coefficients 

1  X  -.2ff,m  N  1  N 

Vm  =  ~fj     2-    ^e    JV    '     for    m  =  ~ "2  +  l  :  2  " 

Thus,  the  DFT  can  be  considered  to  approximate  the  Fourier  transform  t>(wm)  by 
y(wm)  «  AV^. 

We  now  give  a  formal  definition  of  the  DFT1. 


'There  are  several  definitions  of  the  DFT,  as  indicated  in  [Ref.   10]. 
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Definition  A.l  Let  N  be  an  even  positive  integer  and  {vj}  a  sequence2  of  N  complex 
numbers  where  j  =  —  ^  +  1  :  y .  The  discrete  Fourier  transform  of  the  sequence  {vj} 
is  another  sequence  of  N  complex  numbers,  {V^},  whose  elements  are  given  by 

K 

1  2 

Vm  =  ^     £     ^e-V",  (A.l) 

for  in  =  -—  +  1  :  — . 

The  I) FT  may  be  defined  for  N  odd  as  well,  however  we  will  not  need  this  definition 
for  our  discussion.  The  choice  of  indices,  —  y  +  1  :  —  as  opposed  to  0  :  N  —  1,  is 
made  deliberately  in  order  to  simplify  some  of  the  analysis  that  is  to  follow  from 
applying  the  DFT,  and  because  it  is  more  natural,  if  less  familiar.  While  the  choice 
of  indices  does  not  affect  the  applicability  of  the  transform,  care  must  be  taken  to 
avoid  confusion  over  which  grid  points  correspond  to  what  type  of  frequency  (see 
Figure  34).  We  use  the  operator  notation  V{vj]  to  denote  the  DFT  of  the  sequence 
{vj},  and  T>{vj}m  to  indicate  the  mth  element  of  the  transform,  T>{vj}m  =  Vm. 

In  general,  the  output  of  the  DFT,  {Kn},  is  a  complex-valued  sequence.  The 
interpretation  of  the  DFT  coefficients  is  essentially  the  same  as  that  given  for  the 
(continuous)  Fourier  transform.  Also,  Vm  is  even  more  closely  related  to  cm,  the 
Fou tier  series  coefficient.  In  many  ways,  the  DFT  is  more  naturally  viewed  as  an 
approximation  to  cm  than  to  v(u>).  The  rath  DFT  coefficient  Vm  gives  the  "amount" 
of  the  rath  mode  (with  a  frequency  u>m)  that  is  present  in  the  input  sequence  Vj.  In 
contrast  to  the  use  of  modes  of  all  frequencies,  as  in  the  Fourier  transform,  an  ./V-point 
DFT  uses  only  TV  distinct  modes,  with  roughly  y  different  frequencies.  The  modes 
ran  be  labeled  by  the  frequency  index  ra,  and  each  mode  has  a  value  at  each  grid 
point  Xj  where  j  =  — y  +  1  :  y.  Therefore  we  denote  the  jth  component  of  the  rath 


DFT  mode  as 


for  j,  m  =  —  y  +  1 


-jm  -.2ir;m 


-The  terms  sequence  and  vector  are  used  interchangeably  to  denote  the  input/output  of  the  DFT. 
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Figure  34.  Using  centered  indices  (bottom  figure),  the  low  frequencies  correspond  to 
| /.'|  <  j  for  a  single  set  of  sample  points;  the  high  frequencies  correspond  to  |A;|  >  — . 
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Wit  h  non-centered  indices,  low  frequencies  correspond  to  Q  <  k  <  —  and  —  <  k  <  N 
for  a  single  set;  high  frequencies  correspond  to  ^  <  A;  <  ^-. 


The  DFT  allows  us  to  transform  from  the  temporal  (spatial)  domain  into  the 
frequency  domain.  The  Inverse  Discrete  Fourier  Transform.  IDFT,  permits  us 
to  transform  from  the  frequency  domain  back  into  the  temporal  (spatial)  domain.  We 
now  give  a  definition  for  the  IDFT. 

Definition  A. 2  Let  N  be  an  even  positive  integer  and  {Vm}  a  sequence  oj  A  complex 
numbers  where  m  =  —  y  +  1  :  y.  The  inverse  discrete  Fourier  transform  oj  {\m}  is 
another  sequence  of  N  complex  numbers,  whose  elements  are  given  by 


i2ir)m 


(A.2) 


m=-f+l 


for  j  =  --  +  1  :  - 


As  with  the  DFT,  the  IDFT  is  defined  for  /V  odd,  but  we  will  not  need  this  defini- 
tion. Additionally,  we  use  the  operator  notation  V~x  {Vm}  to  denote  the  IDFT  of  the 
sequence  {Vm},  and  V~l{Vm}j  to  indicate  the  jth  element  of  the  inverse  transform. 
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The  discrete  Fourier  transform  and  its  inverse  have  many  useful  properties. 
Below  are  listed  some  of  the  more  important  ones  that  we  will  be  using.  First, 
however,  we  give  the  discrete  orthogonality  property  for  the  complex  exponential, 
which  is  essential  in  working  with  DFTs. 

Property  A.l   Orthogonality.     Let  I  and  m  be  integers  and  let  N   be  a  positive 
integer.    Then 

7V-1 

« — -         i2tv}1  i27rim  * 

£  eV-e — ih  =  JV5N{1  -  m),  (A.3) 

inhere  6]v(k)  '•s  the  modular  Kronecker  delta,  defined  by 

-  J    1     if  k  =  0  or  a  multiple  of  N 

0  otherwise 

Although  we  will  not  specifically  need  it,  we  include  the  relationship  between 
the  DFT  and  the  IDFT: 

Property  A. 2  Inversion.    Let  {vj}   be  a  sequence  of  N  complex  numbers  and  let 
p{c;}  =  {Vm}  be  the  DFT  of  this  sequence.    Then  T>~1  {T>{vj}m}j  —  v3. 

The  remaining  properties  that  we  will  need  are  presented  below.  Proofs  and 
detailed  explanations  are  found  in  [Ref.   10]. 

Property  A.3  Periodicity.    The  complex  sequences  {vj}  and  {Vm}  defined  by  the 
N -point  DFT  pair  (A.l)  and  (A. 2)  are  N -periodic.   That  is, 

Vj  =  vj±n     and     Vm  =  Vm±N     for  all  intergers  j  and  m. 

Property  A. 4  Linearity.     //  {vj}   and  {wj}  are  two  complex-valued  sequences  of 
length  N  and  a  and  j3  are  two  complex  numbers,  then 

V{avj  +  pw3}m  =  aV{vj}m  +  ^{w^m. 

Property  A. 5   Modulation.    If  the  elements  of  the  input  sequence  {vj}  are  multi- 
plied by  ujjn  for  k  a  fixed  integer,  then  the  DFT  of  the  modulated  sequence  is  given 

by 

V{vjUJN}m  =  V{vj}m.k  =  Vm-k- 
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Property  A. 6  Shifting.   If  tin  input  sequence  {vj}  is  shifted  (or  translated)  k  units 
to  ll/t   right,  th(  n 

We  n lake  extensive  vise  oi  the  DFT  in  analyzing  our  solution  process.  The  DFT 
is  an  extremely  useful  tool  in  predicting  the  performance  of  relaxal  ion  schemes  ( ( Chap- 
ter IV)  and  in  evaluating  various  elements  in  the  multilevel  solution  (Chapter  VI)  via 
local  mode  analysis.  Of  the  properties  outlined  above,  the  orthogonality  and  shifting 
properties  are  the  two  that  figure  most  prominently  in  local  mode  analysis. 
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APPENDIX  B.  INTERPOLATION  TOOLS 

Below  are  the  interpolation  tools  developed  by  David  Canright  [Ref.  3]  for 
computing  the  FVE  stencils  using  Maple.  The  indexing  in  the  following  program 
does  NOT  follow  the  indexing  in  the  text:  the  subscript  i  is  the  column  indicator 
and  corresponds  to  the  column  indicator  j  used  in  the  text;  the  subscript  j  is  the 
row  indicator  and  corresponds  to  the  row  indicator  k  in  the  text.  Thus,  the  indices 
( /.  /)  in  this  Maple  program  indicate  the  jth  row  and  the  ith  column,  in  the  text  the 
indices  (h.j)  indicate  the  A;th  row  and  jth  column. 

The  routine  "plane"  returns  the  function  for  the  plane  through  the  three  given 
points.  There  is  no  error  checking,  so  it  has  problems  with  special  cases:  all  3  points 
collinear  (infinite  solutions);  and  vertical  planes  (no  solutions). 

plane  :=  proc(xl,yl,zl,x2,y2,z2,x3,y.3,z3)  local  a,b,c,x,y; 
unapply( 
subs( 

solve(  { 

a  *  x  1  +  6  *  yl  +  c  —  zl, 
a  *  x'2  +  b  *  t/2  +  c  =  z2, 
a  *  x'.l  +  b  *  yS  +  c  =  z'S  }, 
{a,6,c}), 
a  *  x  +  b  *  y  +  c), 

x,y  ); 

end: 

To  interpolate  the  indexed  function  z  on  the  6  triangles  surrounding  the  point 
(/'  *  /*,  j  *  h),  use  the  6  functions  (counter  clockwise  from  northeast)  "tri:  ene,  nne, 
tiw,  wsvv,  ssw,  se".  Note:  these  return  functions  of  two  variables  (for  example  r.z). 
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triene  :=  (ij,z)  -  >  plane(i*h,j*h,z._p,(i+l)*h,j*h,z._e,(i+l)*h, 

(j  +  l)*h,z._ne): 
triune  :=  (i,j,z)  -  >  plane(i*h,j*h,z._p,i*h,(j+l)*h,z._n,(i+l)*li, 

(j+l)*h,z._ne): 
trinw  :=  (i,j,z)  -  >  plane(i*h,j*h,z._p,i*li,(j+l)*h,z._n,(i-l)*h, 

j*h,z._w): 
triwsw  :  =  (i,j,z)  —  >  plane(i*h,j*h,z._p,(i— l)*h,j*h,z._w,(i—  l)*h, 

(j-l)*h,z._sw): 
trissw  :-  (i,j,z)  -  >  plane(i*h,j*h,z._p,i*h,(j  —  l)*h,z._s,(i-l)*h, 

(j-l)*h,z._sw): 
trise  :=  (ij,z)  -  >  planef i*hj*h.z._p,i*h,(j  — l)*h,z.js,(i+l)*h, 

j*h,z._e): 


To  convert  the  notation  from  subscripted  back  to  indexed,  use  utoindex(expression, 
variable)",  where  the  variable  is  to  become  indexed. 

toindex  :=  proc  (expr,var)  local  n,m; 

subs( 

zip(  (x,y)  -  >  (x=y), 

[var.  (_sw,_s,_se,_w,_p,_e,_nw,_n,_ne)], 

[seq  (seq  (seq  (var[i+n,j+m],  n=  —  l..l),m=  — 1..1)]), 

expr  ); 

end: 


Now  we  define  shorthand  for  limits  of  integration  about  the  general  point  (i,j). 
Here,  rdiag  means  the  diagonal  line,  where  r  depends  on  z.  (This  assumes  integrating 
first  in  7-,  then  in  z.) 

rp  :=  i*h;  rw  :=  (i-l/2)*h;  re  :=  (i+l/2)*h;  rdiag  :=  i*h-(z-j*h); 
zp  :=  j*h;  zs  :=  (j-l/2)*h;  zn  :=  (j+l/2)*h; 

Now  compute  the  volume  integral  for  the  unknown  temperature  T  (which  is  a 
surface  integral  of  T  after  factoring  out  the  27T  from  the  ip  integration): 


int  (int  (triene  (i,j,T)(r,z)*r,  r=rdiag..re),z=zp..zn)  + 
int  (int  (trinne  (i,j,T)(r,z)*r,  r=rp..rdiag),z=zp..zn)  + 
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int  (int  (trinw  (i,j,T)(r,z)*r,  r=rw..rp),z=zp..zn)  + 
int.  (int  (triwsw  (i,j,T)(r,z)*r,  r=rw..rdiag),z=zs..zp)  + 
int  (int.  (trissw  (i,j,T)(r,z)*r,  r=rdiag..rp),z=zs..zp)  + 
int  (int  (trise  (i,j,T)(r,z)*r,  r=rp..re),z=zs..zp); 

factor!  simplify(" )); 
toindex(",T); 


This  set  of  computations  produces  the  stencil  entries  for  the  volume  integral: 

—  ((32»  +  UjTi+ij  +  (32.  +  5)riti_,  +  224xi;-j  +  (32t  -  1 1  )Tt.hj 

.,s  i 

(16*  -  o)Ti-U-i  +  (32.  -  5)TM+1  +  (16t  +  5)TI+l,J+1) . 

In  similar  fashion,  the  surface  integral  of  VT  is  computed  (which  is  a  circula- 
tion integral  of  T  after  factoring  out  the  27r  from  the  tp  integration): 

int  (  +  D[1]  (triene(i,j,T))(re,z)*re,z=zp..zn)  + 
int  (  +  D[2]  (trinne(ij,T))(r,zn)*r,r=rp..re)  4- 
int  (  +  D[2]  (trinw(i,j,T))(r,zn)*r,r=rw..rp)  + 
int  (-D[l]  (trinw(i,j,T))(rw,z)*rw,r=zp..zn)  + 
int  (-D[l]  (triwsw(i,j,T))(rw,z)*rw,z=zs..zp)  + 
int  (-D[2]  (trissw(i,j,T))(r,zs)*r,r=rw..rp)  + 
int  (-D[2]  (trise(i,j,T))(r,zs)*r,r=rp..re)  + 
int  (  +  D[1]  (trise(i,j,T))(re,z)*re,z=zs..zp); 

factor(simplify(" )); 
toindex(",T); 


This  set  of  computations  produces  the  stencil  entries  for  the  surface  integral: 


107 


108 


REFERENCES 


Stephen  F.  McCormick.  Multilevel  Adaptiv(  Methods  for  Par/ml  Differential 
Equations.  Society  for  Industrial  and  Applied  Mathematics,  Philadelphia,  1989. 

Achi  Brandt.  Multigrid  Ticliiuqtu  s:  l!)S.j  (luidt  with  Applications  to  Fluid  Dy- 
namics. The  von-Karman  Institute  for  Fluid  Dynamics.  Rhode-Saint-Genese, 
Belgium.  1984. 

[3]  William  L.  Briggs.  .4  Multigrid  Tutorial.  Society  for  Industrial  and  Applied 
Mathematics,  Philadelphia,  1987. 

David  Canright.  Thermocapillary  effects  on  weld  pool  shape,  working  draft. 
unpublished,  1994. 

II.  S.  Carslaw  and  J.  C.  Jaeger.  Conduction  of  Heat  in  Solids,  Second  Edition. 
Oxford  University  Press,  Oxford,  Great  Britain,  1959. 

Eugene  Isaacson  and  Herbert  Bishop  Keller.  Analysis  of  Numerical  Methods. 
Dover  Publications,  Mineola,  New  York,  1994. 

G.  D.  Smith.  Numerical  Solution  of  Partial  Differential  Equations:  Finite  Dif- 
ference Methods.  Oxford  University  Press,  Oxford,  Great  Britain,  1985. 

David  Canright.  Interpolation  tools,  for  construction  of  stencils  for  FVE  method, 
personal  communication,  1994. 

Gene  H.  Golub  and  Charles  F.  Van  Loan.  Matrix  Computations.  The  .Johns 
Hopkins  University  Press,  Baltimore,  1983. 

10]  William  L.  Briggs  and  Van  Emden  Henson.  The  DFT:  an  Owners  Manual  for 
the  Discrete  Fourier  Transform.  Society  for  Industrial  and  Applied  Mathematics, 
Philadelphia,  to  appear. 

11]  C.  Liu.  The  Finite  Volume  Element  (FVE)  and  Fast  Adaptive  Composite  Grid 
Methods  (FAC)  for  the  Incompressible  Navier-Stokes  Equations.  Copper  Moun- 
tain Conference  on  Multigrid  Methods,  1(1):  1—21 ,  1993. 

12]  K.  W.  Hamming.  Numerical  Methods  for  Scientists  and  Engineers,  Second  Edi- 
tion. Dover  Publications,  Mineola,  New  York,  1986. 

[13]  David  Canright.  Thermocapillary  Flow  Near  a  Cold  Wall.  Phys.  Fluids, 
6(4):1415-1424,  1994. 


109 


110 


INITIAL  DISTRIBUTION  LIST 


1.  Defense   Technical  Information  Center 
1  ameron  Station 

Alexandria,  Virginia  22304-6145 

2.  Library,  ( 'ode  52 

Naval  Postgraduate  School 
Monterey.  California  93943-5101 

3.  Professor  Richard  Franke,  Code  MA/Fe 
Department  of  Mathematics 

Naval  Postgraduate  School 
Monterey,  California  93943-5216 

I.    Professor  David  Canright,  Code  MA/Ca 
Department  of  Mathematics 
Naval  Postgraduate  School 
Monterey,  California  93943-5216 

5.  Professor  Van  Emden  Henson,  Code  MA/Hv 
Department  of  Mathematics 

Naval  Postgraduate  School 
Monterey,  California  93943-5216 

6.  Professor  Donald  Danielson,  Code  MA/Dd 
Department  of  Mathematics 

Naval  Postgraduate  School 
Monterey,  California  93943-5216 

7.  Professor  Christopher  Frenzen,  Code  MA/Fr 
Department  of  Mathematics 

Naval  Postgraduate  School 
Monterey,  California  93943-5216 

(S.    Director,  Training  and  Education 
MCCDC,  Code  C46 
1019  Elliot  Road 
Quantico,  Virginia  22134-5027 

9.   Requirements  Division 
MCCDC,  Code  C442 
2042  Broadway  Street,  Suite  11 
Quantico,  Virginia  22134-5021 


111 


1  .i^ARY 

•  KADUATE  SCHOOL 
MONTEREY  CA   83043-5101 


SAYIOPD  s 


